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This  document  supersedes  the  previous  AFOSR  Scientific  Report  cited 
below.  The  present  report  contains  a  number  of  improvements  in  the  numerical 
methods  used,  and  is  applicable  to  a  wider  range  of  blade  geometries.  In  addi¬ 
tion,  a  number  of  errors  in  Ref.  1  have  been  corrected,  and  more  complete 
descriptions  are  given  for  certain  features  of  the  method. 

The  author  is  very  indebted  to  Joseph  P.  Nenni,  John  R.  Moselle,  and 
Marcia  J.  Williams  for  their  assistance  in  the  development  of  this  program. 
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ABSTRACT 


This  report  contains  the  description  of  a  computer  program  for  evalu¬ 
ating  the  Ives  transformation,  which  maps  a  cascade  of  turbine  or  compressor 
blades  conformally  into  a  rectangle. 
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Section  1 
INTRODUCTION 


Some  of  the  most  important  recent  advances  in  computational  fluid 
dynamics  have  been  made  possible  by  the  use  of  algorithms  which  solve  the 
equations  of  motion  in  boundary-conforming  coordinates.  Thus  the  development 
of  methods  for  carrying  out  these  coordinate  transformations  has  received 
considerable  attention  (see,  for  example.  Reference  2).  The  methods  in  current 
use  can  be  divided  into  three  categories:  those  which  apply  algebraic  shearing 
or  stretching  of  the  coordinates,  those  which  are  generated  numerically  by 
solving  a  Poisson  equation,  and  those  which  are  based  on  a  conformal  transforma¬ 
tion. 


This  report  is  concerned  with  one  example  from  the  third  category, 

3 

namely  the  transformation  introduced  by  Ives  and  Liutermoza.  It  is  capable 
of  transforming  a  cascade  of  airfoils  into  a  rectangle,  thus  facilitating  the 
application  of  the  solution  algorithms  mentioned  above.  The  content  of  this 
report  is  a  review  of  the  transformation  itself,  together  with  a  practical 
numerical  procedure  for  applying  it  to  a  given  cascade.  This  procedure  in¬ 
volves  a  number  of  choices,  involving  branch-point  locations,  tolerances  for 
various  iteration  sequences,  and  formulas  for  the  evaluation  of  certain  special 
functions. 


Section  2  below  contains  a  description  of  the  transformation  itself, 
including  modifications  that  enable  the  mapping  of  blades  having  a  rounded 
trailing  edge,  aid  the  calculation  of  grids  in  which  one  of  the  coordinate  lines 
joins  the  trailing  edge  to  the  point  at  downstream  infinity.* 


2.  Thompson,  J.F.,  ed..  Numerical  Grid  Generation,  Elsevier  Science 
Publishing  Co.,  New  York  (1982). 

3.  Ives,  D.C.  and  Liutermoza,  J.F.,  "Analysis  of  Transonic  Cascade  Flow 

Using  Conformal  Mapping  and  Relaxation  Techniques",  AIAA  Journal  15 
(1977),  pp.  647-652.  — 

4.  Rae,  W.J.,  "Modifications  of  the  Ives- Liutermoza  Conformal -Mapping 
Procedure  for  Turbomachinery  Cascades",  ASME  Paper  83-GT-116  (March  1983). 
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Section  2 

CONFORMAL  TRANSFORMATION  METHOD 


The  method  of  Ives  §  Liutermoza  consists  of  a  sequence  of  transforma¬ 
tions,  which  map  a  two-dimensional  cascade  into  a  rectangle.  The  notation  and 
coordinate  system  used  to  define  the  cascade  are  shown  in  Fig.  1.  The  %  -co- 


The  quantities  3  ,  n.  and  the  angle  Y  denote  the  "streamwise,  normal"  coordinates, 

in  terms  of  which  the  blade  profiles  are  sometimes  defined.  These  reduce  to 

* 

the  tC  ,  ^  set  if  H  is  taken  as  zero.  The  origins  of  both  of  these  coordinate 
systems  are  arbitrary. 

These  coordinates  define  a  complex  variable  2  : 

2  -  S  +  in, 

The  points  ZN  and  ZT  are  taken  anywhere  near  the  centers  of  curvature  of  the 
leading  and  trailing  edges,  while  ZLE  and  ZTE  are  points  which  divide  the 
"pressure"  side  of  the  blade  (i.e.,  its  concave  surface)  from  the  "suction" 
side  (its  convex  surface) .  These  points  can  be  chosen  anywhere  on  the 


*  .  +  _9 

Actually,  n  must  not  be  set  identically  equal  to  zero,  but  to  the  value  -  10  , 

where  the  ±  signs  apply  to  compressors  or  turbines,  respectively. 
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leading-  and  trailing-edge  contours;  ZTE  is  the  point  that  will  be  connected 
to  the  "point"  at  downstream  infinity  by  one  of  the  grid  lines,  if  that  option 
is  chosen  (i.e.,  ISHEAR=1). 

For  the  case  of  a  sharp  trailing  edge  (ITE=0),  ZT  must  equal  ZTE, 
and  for  a  sharp  leading  edge  (ILE=0),  ZN  must  equal  ZLE.  The  included  angle 
at  a  sharp  trailing  edge,  X  ,  must  be  specified.  This  is  illustrated  in 
Fig.  2,  for  the  cascade  used  by  Rae  and  Homicz.^ 


Figure  2.  Blade  Geometry  of  Ref.  5 

This  blade  row  has  a  slant- gap/ axial  chord  ratio  SG/Ca»  i  ,  where  SCr  is  the 
slant  gap: 


5.  Rae,  W.J.,  and  Homicz,  G.F.,  "A  Rectangular-Coordinate  Method  for 
Calculating  Nonlinear  Transonic  Potential  Flowfields  in  Compressor 
Cascades",  AIAA  Paper  78-248  (January  1978). 
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L 


SG  = 


(2-2) 


and  a  stagger  angle  i  of  32.91°.  This  geometry  is  used,  below,  to  illustrate 
the  steps  in  the  Ives  transformation.  A  comparable  set  of  illustrations,  for 
a  cascade  of  turbine  blades  with  rounded  trailing  edges,  is  given  in  Ref.  4. 

Figure  2  shows  the  relation  between  the  trailing-edge  angle  X  and  an 
exponent  EX  (which  is  used  in  one  of  the  transformation  steps  described  below) . 
This  relation  applies  only  for  a  sharp  trailing  edge.  For  a  rounded  trailing 
edge,  Ex  is  not  related  to  the  180-degree  trailing-edge  angle,  but  is  chosen 
as  a  number  in  the  range  0.2  to  0.4,  as  described  below. 

The  blade  shape  is  input  as  two  tables  of  coordinate  pairs,  one  for 
the  pressure  surface,  and  one  for  the  suction  surface.  These  coordinates,  plus 
the  leading-  and  trailing-edge  points  2LE  and  ZTE  are  then  arranged  in  an  array 
indexed  by  KJ,  where  KJ=1  at  the  trailing  edge,  and  where  the  numbering  pro¬ 
ceeds  around  the  pressure  side  to  the  leading  edge  (KJLE) ,  and  then  along  the 
suction  side  to  the  trailing  edge,  where  the  point  denoted  by  KJMX  is  a  repeat 
of  that  denoted  by  KJ=1.  This  notation  is  shown  in  Fig.  3. 


*  j-  m 


ZLE  - 5 

(KJ  =  KJLE) 


ZTE 


5  ( KJ  *  1  s!  KJMX) 


Figure  3.  Notation  for  Blade-Surface  Coordinates 

The  quantities  KJS  and  KJP  need  not  be  equal;  they  are  limited  to  a  maximum 
value  of  80  by  a  dimension  statement  in  the  current  version  of  the  program. 


The  first  step  in  the  Ives  transformation  is 
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~  I*  Hill 

1  «  +  CG  ] 

•  f_  Z-2.  ‘ 

6-VL  4  77  — - — —  ^ 

I  + 


The  fact  that  only  differences  of  2  -values  are  used  is  what  account  or  the 

arbitrariness  of  the  origin  in  Fig.  1.  On  the  suction  (  5  )  and  pre  -e  (  P  ) 

sides,  the  function  ^(2)  has  the  Fortran  equivalents:: 


RDS  ( K )  e 
R  DP  CK)  e 


i  rus  ( K ) 

L  THP(K) 


(2-4) 


as  follows: 


The  arguments  of  the  sine  functions  can  be  written  in  a  simpler  form, 
i:  „  AZ(//*(5j 


*  .  r  * 


H 


*  — —  \fU  A2  •  t  y  +  Jl**s  4z  •  Cm  i 
SQ  l 

+  t  £  d-rrt  A  2  •  >^v>  y  —  Rx  £  a  ' »']} 


where  At  stands  for  2 ~2^  or  2“2„  in  the  expressions  for  and  <i^3  ,  respectively. 
By  noting  that  (see  Fig.  4) 


X  -  X 


«  5  c*o  X  ~  n  AA+x  y 


f  ~  y  +  n  cm  y 


Figure  4.  Coordinate  Relations 
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it  follows  that 


=  j^r  |(S  -  ST)  Ac**  i  +  (n -nT)  c+o  i  +  c  -n.r)  -  (s  ~sr)  c**g 


(2-6 


Similarly, 


A  =  —  j 

3  SCr  L 


-  L  (X 


(2-7; 


The  sine  transformations: 


■nap  the  strip  -  tt/2  $  4  +  71/2  into  the  entire  (or  ^  ) 

plane,  with  cuts  along  the  real  axis  from  -«»  to  -1  and  from  +1  to+<»; 


Figure  S.  The  Mapping  jr 


In  order  to  enforce  this  single- valuedness,  the  argument  of  the  sine  function 
should  be  displaced  by  the  appropriate  multiple  of  H  +•  i  &  ,  whenever  the 
real  part  of  its  argument  falls  outside  the  strip  noted  above,  i.e.,  if 
\RtZ,y\  >  -j-  ,  then  set  £  »  (Z,  +  nTT)-nTr,  where  n.  is  chosen  so  as  to  make 
I  Re  ( %  *  nn  )  I  c  -j- .  It  turns  out,  however,  that  this  displacement  of  the 
argument  need  not  be  made  explicitly,  when  evaluating  the  complex  sine.  This 
can  be  seen  by  noting  that 

.  -  to  ,  . 

(  Re  )  *  (  R  C^O-4  d  )  (  R  sQjw  0 )  +  C  CO<  ( R  C*a  9  )  (  R  ) 

and  that 
,  ^  £ 

[(Re  +  n  Tt)  -  n  TC  ]  m  Oo  r>  rt  •  ( R  Q  +  n  It  -4-  L  R  &) 

a  C  R  C**  0  )  ( R  AM"  6  )  +  L  (  R  C*0  9  )  (  R  A4*r>  Q) 


Thus,  no  special  treatment  is  required  for  the  argument  of  the  complex  sine. 
The  result  of  this  transformation  is  shown  in  Fig.  6,  where  S  and  P  denote 
the  suction  and  pressure  sides.  The  quantity  ej(£)  is  then  formed  from  the 
ratio  of  the  sines;  this  curve  is  shown  in  Fig.  7.  The  interior  of  the  blade 
lies  outside  the  closed  curve  in  this  plane.  Because  this  function  is  to  be 
raised  to  a  power,  its  argument  must  be  defined  on  the  suction  and  pressure 
sides  in  such  a  way  that  the  exponentiation  will  map  the  trailing-edge  region 
into  a  straight  segment.  In  the  computer  program,  these  arguments  are  first 
found  from  the  Fortran  DATAN2  function,  which  returns  angles  in  the  range  from 
-71  to +7 1  .  These  arguments  are  then  adjusted,  by  defining  the  argument  for 
KJ=2,  and  by  then  examining  the  points  KJ=3  to  KJMX,  adding  or  subtracting 
2  TV  to  the  argument  whenever  the  cut  (along  the  negative  real  axis)  in  the 
DATAN2  function  is  traversed  counterclockwise  or  clockwise,  respectively.  The 
arguments  defined  for  KJ=2  use  the  following  convention:  for  H  >0  (i.e. ,  a 
stagger  corresponding  to  a  compressor  blade  row)  the  argument  at  KJ=2  is  taken 
to  be  negative.  For  all  other  cases,  the  argument  at  KJ=1  is  accepted  as 
returned  by  the  DATAN2  function. 
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"The  Transformations  and  ^ 


The  next  step  is  the  exponentiation,  defined  as: 

yK 

ci  -  [gC2)l  :  K  •  2  -  —  ~  exinv  h  —  C2-8) 

*  rr  £X 

The  result  of  this  step  is  shown  in  Fig.  8.  As  noted  above,  the  value  of  ffx 
for  a  sharp- trailing -edge  blade  is  fixed  by  the  angle  X  ,  whereas  for  a  round 
trailing  edge  a  range  of  values  of  £X  ,  in  the  neighborhood  of  0.3,  will  produce 
a  curve  in  the IL  -plane  that  is  "star-shaped",  i.e.,  its  polar  coordinates  will 
have  a  radius  that  is  a  single-valued  function  of  the  angular  coordinate. 


The  next  transformation  is  a  bilinear  one,  whose  purpose  is  to  produce 
a  curve  in  the  to -plane  that  can  be  mapped  into  a  unit  circle  in  a  subsequent 
step: 


oj-a. 
OJ  -b 


CL -b 


A. 


A 


Ui 


1  - 


-CL 

C 


(2-9) 


where  a- ,  b  ,  and  C  are  complex  constants.  Three  conditions  must  be  assigned, 
in  order  to  evaluate  these  constants:  Ives  suggests,  for  two  of  them,  that 
the  images  of  upstream  and  downstream  infinity  be  mapped  into  u 1  .  As 
the  third  condition,  he  recommends  that  the  centroid  of  the  blade- surface  image 
in  the  cO -plane  be  forced  to  be  close  to  the  origin.  This  condition  was  applied 
in  Ref.  1;  however,  it  has  been  found  simpler  to  impose  a  condition  on  the 
ratio  between  the  maximum  and  minimum  radii  in  the  uJ -plane,  as  outlined  below. 
The  calculation  of  the  centroidal  location  has  been  retained  in  the  present 
code,  for  informational  purposes.  (Details  on  how  this  is  calculated  can  be 
found  in  Ref.  1). 


The  locations  of  the  images  of  upstream  and  downstream  infinity  in 
the  fl  -  andil  -planes  can  be  expressed  as  follows:  in  general. 


fl.2) 


AA-y'.  (  fit-  %t)  CAaSi XJUr>  ) 

t  at  - 


ii 


The  Mapping  J1CZ) 


As  Z 


-  00 


along  line  ^  of  Fig.  9, 


m  finite j  yjm  -+  +  M  *  0  ,  Jw  £}  — *  +  m 

Using  the  large-argument  approximations  for  the  hyperbolic  functions  then 
gives 

-«■  <&-?,  x*  ***•* 

r-°°;  -  *  c  *  c  c  (2-10) 


where  a  little  algebra  reveals  that 

%  - /<  *  M,(z;i  - z;3)  *  —  {a^2t  -**)  &W-  Jvw  rzr- zw;  ^v.ir'} 

x*-*  s  &.<*,;  *  -  +  Jm(z7-z*; 

As  H +  "  along  line  of  Fig.  9, 

%1  *  ^  >  **  m  h'nite  ,  >J*» •  2J}  — *  -«» 

and  these  lead  to 


(4-  eo) 


tsf  -?,)  -  <•  A(^s; 


where  £Oj,)  =  -  4L(S,) 


and 

Uk  $  1%  4  »/k 

il’  *  [^  <- <*>]  =  «  e  ,  /I  =  [<j  (+«»)] 

Note  that 

Jl*  fl”  *  ; 


-x 

e  e 


(2-12) 
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The  constants  Ol  and  b  can  be  expressed  in  terms  of  C  by  the  two  equations: 

-1-a-  .  - 

C  -  =  A 

-1-b 


The  solution  is: 


1-0-  _  * 

C  -  2  SL 

1-b 


-  2  C  +  SL~  +  SL+  ,  ,  ,  SL 

- - - - -  ;  a.  »  1  -  (1  -  b  )  — 

SL-  SL  C 


-  Ec  +  F 


where  F  and  F  are  known  quantities: 


si*  -nr 


,  F  2 


(2-13) 


-a*  -id" 
si*  -  XI" 


(2-14) 


In  terms  of  the  parameters  E  ,  F  ,  and  c  ,  the  transformation  can  be  written 


SI  = 


E  -  Fo  +  (Ec  -F)  SI. 
C  -  SL 

(ui*F)c  -  E 
oo-F  +  E  c 


(2-15) 


(2-16) 


The  latter  relations  can  be  used  in  an  iteration  procedure  to  find  a  value  of 
C  that  will  minimize  the  ratio  RMAX/RMIN,  where  RMAX  and  RMIN  are  the  maximum 
and  minimum  values  of  |cd I  over  the  set  defined  by  KJ=1  to  KJMX.  The  iteration 
process  is  as  follows:  on  alternate  iterations,  values  of  C  are  chosen  that 
will  either  reduce  RMAX  or  increase  RMIN.  This  is  done  by  solving  Equation  2-16 
for  C  : 

SLcJ  +  B  -  F  SL 


uJ  +  F  -  E  SL 


(2-17) 


On  the  first  iteration,  C  *  f  +  4.  is  used;  the  values  of  RMIN  and  the  index 
KJ=KJMN  at  which  it  occurs  are  then  found.  Next,  a  new  value  of  C  is  calculated. 


15 


using  Eq.  2-17,  such  that  the  new  value  of  0J  (KJMN)  will  equal  1.1  times  the 
value  just  found.  For  this  new  mapping  into  the  60-plane,  the  value  of  RMAX 
and  the  index  KJ=KJMXX  at  which  it  occurs  are  found.  For  the  third  iteration, 
a  value  of  C  is  used  such  as  to  generate  a  new  value  of  cO  (KJMXX)  equal  to 
0.9  times  the  previous  one.  This  alternating  cycle  is  then  continued  until 
the  ratio  RMAX/RMIN  is  less  than  the  tolerance  RTOL.  This  tolerance  is  assigned 
a  default  value  of  3.0;  values  up  to  6.0  have  been  handled  successfully  by  the 
subsequent  steps  in  the  transformation. 

The  iteration  on  C  can  be  bypassed,  if  C  is  already  known,  by  setting 
IG0T=1  and  reading  in  the  value  of  C  .  Figure  10,  reproduced  from  Ref.  1,  was 
generated  this  way. 


Next,  the  blade-surface  image  in  the  u/-plane  is  mapped  into  the  unit 
circle  in  the  c,  -  plane  with  the  trailing  edge  at  £,  *  1  by  a  variant  of  the 
Theodorsen- Garrick  transformation: 

uj  =  %  jupp  i  £  +  1  * 

l  ^*o  * 

To  determine  the  coefficients,  the  values  of  to  and  (,  on  the  blade  surface 
are  written  as 

19  „ 

60  =  r  (.9)  C  ,  £  =  1  tf  (2-19) 


} 


(2-18) 


Then  the  real  and  imaginary  parts  of  the  transformation  are 


Jk> i  r  * 

Q  - 


*o  +  2  ^  j  4*  ~  i 

i*' 

+  [flj  ^  }  $] 


(2-20) 


(2-21) 
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The  coefficients  are  then  determined  by  the  following  iteration  procedure: 
an  equally- spaced  array  of  values  of  ^  is  set  up,  and  all  the  coefficients 
are  initially  set  equal  to  zero.  Then  the  second  equation  gives  &  =  <p  as 
the  first  approximation  for  6  .  For  each  of  these  values  of  Q  ,  a  corres¬ 
ponding  value  of  f  can  be  found,  from  a  spline  fit  to  the  coordinates  of 
the  blade-surface  image  in  theu>-plane.  These  known  values  of  i»i("  are  then 
used  in  the  first  of  the  equations  above,  to  find  the  next  approximation  to 
the  and  8 ^  coefficients .  These  coefficients  can  then  be  used  to  give  the 

next  approximation  to  0(<fr) ,  and  the  process  is  continued  until  convergence 
is  reached,  to  some  preassigned  tolerance.  Fast  Fourier  transform  tech- 
techniques6  can  be  used  in  the  processes  of  evaluating  the  second  equation,  for 
known  values  of  the  coefficients,  and  of  determining  the  coefficients  from 
the  first  equation  with  known  values  of  Jbn  r  .  These  techniques  were  applied 

in  the  present  calculations.  The  details  are  given  in  Appendix  A. 


Sufficient  conditions  for  convergence  of  this  iteration  process 
have  been  discussed  by  Warschawski. 7  For  the  present  case,  these  conditions 
are  not  met;  in  particular,  it  is  required  that  the  maximum  and  minimum  values 
of  r  (6)  obey  the  relation: 


\ 


RMflX ' 


RmiN 


RMfiX 

<  0.295  or  -  <  7.4>7 3 

RMIN 


This  condition  is  not  met  in  the  present  case,  where  RMAX/RMIN  is  approxi¬ 
mately  2.S.  It  was  found,  however,  that  the  iteration  process  would  converge 


6.  Cooley,  J.W.,  Lewis,  P.A.W.,  and  Welch,  P.D.,  "The  Fast  Fourier  Transform 
Algorithm:  Programming  Considerations  in  the  Calculations  of  Sine, 

Cosine  and  Laplace  Transforms",  Journal  of  Sound  and  Vibration  12, 

(1970)  pp.  315-337 . 

7.  Warschawski,  S.E.,  "On  Theodorsen's  Method  of  Conformal  Mapping  of 

Nearly  Circular  Regions",  Quarterly  of  Applied  Mathematics  3,  (1945) 

pp.  12-28. 
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if  a  relaxation  parameter  was  used,  i.e.,  the  values  of  6  called  for  by  the 
second  equation  [  called  were  not  used  in  the  first  equation,  but  were 
replaced  by  B  new  =  0.10*  +  0.9  Q  ,.  With  this  relaxation  factor,  the 
iterations  were  convergent:  after  68  iterations,  the  maximum  change  in  any 
of  the  values  of  6  was  less  than  10  radians.  The  variation  of  d  with  <p  is 
shown  in  Fig.  11.  Calculations  for  other  cases,  not  shown  here,  have  required 
relaxation  factors  as  low  as  0.02  for  convergnece.  A  recent  review  paper  by 
Henrici  (Reference  8)  calls  attention  to  the  applicability  of  under-relaxation 
in  this  problem. 

The  next  transformation  is 


4-/3 


(2-22) 


where  ct,  /3,  and  are  chosen  so  as  to  place  the  images  of  2  =  ±  eo  at 
=  +  S  ,  while  the  blade  surface  continues  to  be  the  unit  circle .  These  images 

are  located,  respectively,  at  t,  and  £  *  .  Explicit  formulas 

*  3 

for  cc.  ,  /3 ,  and  i  in  terms  of  and  are  given  by  Ives  as 


2  -  1 1  +  2  |  K*  I 
Kfl  -  <sl2 


{Vi  %  +  Vz*-}'  r  ,  Vi  x  -  vi*-  /'i 1 1 

«  £.  ^  ~  As)  ~  ^  *  4  a ,)  3  / 

•«5,a  ^  "  4®  )  +  (4*  +  ~2  / 4* 


8.  Henrici,  P.,  "Fast  Fourier  Methods  in  Computational  Complex  Analysis". 
SIAM  Review  21,  (1979)  pp.  481-527. 
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Figure  11.  Variation  of  &  vs.  $  on  the  Blade 


1 


20 


I 

I 


fl  - 


Zf,  +  i,e  ~  2  oc 


i  • 


<,»  - 


(2-23) 


Mokry  (Reference  9)  has  pointed  out  that  the  formula  for  S can  be  simplified, 
as  follows:  define 


f 

If 


Then 


C  5 


t**  ~  1 

*  4s 


ICI 


-  Vi  ci  *~-Y  ,  i  = 


a  * 


'AICL  ,  & 

C  -  £  ici 


c  -  Sici_  X>B 

I  Cl  Kg'SC 

ICI  - Sc  <3 
ICI  i;e  -  SC 


(2-24) 


The  computer  program  described  below  does  not  use  these  simplifications. 


[ 

I 

I 

n 


Next,  it  is  necessary  to  find  and  4S  ,  given  u>  m  -  1.  This 
was  done  by  Newton- Raphson  iteration: 


<f\4 


In) 


%  - 


QCQ 
o LG 


o U, 


9.  Mokry,  M.,  "Comment  on  Analysis  of  Transonic  Cascade  Flow  Using  Con¬ 
formal  Mapping  and  Relaxation  Techniques",  AIAA  Journal  16,  No.  1, 
(January  1978)  p.  96. 


n 
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L 


where 


In  order  to  find  an  initial  guess  for  5  ,  a  preliminary  calculation  was  made, 
for 

151  =  0.49  (.1)  0.99,  arg£  =  -60°  (10°)  +  60° 

From  this  set,  the  value  of  ^  which  gave  u) nearest  to  +1  was  chosen  as  the 
initial  guess.  This  set  was  then  repeated  with  5  replaced  by  -  %  ,  to  obtain 
the  value  of  ^  that  gave  oi  nearest  to  -1.  The  locations  of  key  points  in  the 

Z,  and  q  planes  are  shown  in  Fig.  12. 

When  the  values  of  Z,n  and  are  known,  the  blade- surface- image  points 

can  be  mapped  from  the  Z,  -plane  to  the  q  -plane.  In  both  planes,  these  points 

are  located  on  unit  circles,  so  it  is  only  necessary  to  interpolate  in  Fig.  12 
to  find  the  Z,  -values  for  the  points  defining  the  blade  surface.  This  is  done 
with  a  call  to  the  spline-fit  subroutine,  and  the  values  returned  by  it  are 

replaced  by  linear  interpolation  in  regions  where  the  &  ,  curve  is  so  steep 

that  the  spline  fit  returns  non-monotonic  values. 

The  final  transformation  now  uses  an  elliptic  function  to  map  the 

unit  circle  in  the  rj  -plane  (with  cuts  along  the  real  axis  from  ±5  to  the 
circle)  into  a  rectangle: 

q  *  3  Sn  (  |  ,4)  (2-25) 

where  the  parameter  k  is  given  by 

•k  -  SZ  (2-26) 
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<  -  PLRtie.  r\  -  Plrk>£ 


Figure  12.  The  Z,  and  Planes 


The  inverse  of  this  transformation  is 


This  latter  transformation  is  used  to  find  the  images,  in 


(2-27) 

/v 

the  ^  -plane,  of 


the  blade-surface  points.  This  requires  an  expression  for  the  real  and  imag¬ 


inary  parts  of  the  incomplete  elliptic  integral  of  the  first  kind;  convenient 
formulae  for  this  purpose  are  given  by  Nielsen  and  Perkins^  who  show  that. 


if 


q  -  S  C  T  +  i  S  ) 


(2-28) 


then 

where 


\  ^cFCVF1  ,  *') 

-k'  -  V/- a1'  -  v? 


(2-29) 

(2-30) 


10.  Nielsen,  J.N.  and  Perkins,  E.W.,  "Charts  for  the  Conical  Part  of  the 
Downwash  Field  of  Swept  Wings  at  Supersonic  Speeds",  NACA  Technical 
Note  1780,  (December  1948),  Appendix  C. 
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and  where  F  (  •  ■  •  )  denotes  the  incomplete  elliptic  integral  of  the  first 
kind,  with  real  arguments  A  and  <r  given  as  functions  of  T ,  6  and  A  : 


A  =■  [/  +Z*+  62-  V(l  -T2)*  +  <?*(<?*+  2  +  2T*/ ]  [;  +  A2  Ct2+  £*) 

-  V(7-  4*rV  1  ±*6\Z  +  Z±2TZ  + 

c r  *  [r2  +  <52- *]/[ra  +  <y4_  x  +  1  -  (**+£*)] 


These  formulae  are  equivalent  to  those  following  Eqn.  115.01  of  Byrd  and 
Friedman;11  the  formula  for  A  has  been  rearranged  slightly  from  the  form 
given  by  Nielsen  and  Perkins,  to  avoid  the  occurence  of  negative  values 
under  the  square  root  sign,  which  can  sometimes  happen  in  the  numerical  evalu¬ 
ations  when  6  =  0  .  These  formulas  are  correct  along  the  branch  cuts, 
where  <5  *  0  and  |  qj  /■S  >  1  . 


Numerical  evaluations  of  these  elliptic  integrals  were  done  using 

12 

twelve  terms  in  the  formulae  of  Luke: 

_  W  <j>  —  t^bnn  U 


?,.(*<*)  -  jr  \*  +  zi f 

c  L  m  - 1 


$  ta^n  ’  (  <Tm  <f>) 


2 


] 


(2-32) 


11.  Byrd,  P.F.,  and  Friedman,  M.D.,  Handbook  of  Elliptic  Integrals  for 
Engineers  and  Physicists.  Springer  Verlag,  Berlin  (1954). 

12.  Luke,  Y.L.,  "Approximations  for  Elliptic  Integrals",  Mathematics  of 

Computation,  22^  (1968)  627-634. 
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where 


(p  <  Tt/2  ,  <TW  =  yi~-kz  a^2  0^,  ,  0  =  nr\TC/25 


For  the  case  where  <p  -  TC/2  (the  complete  integral)  the  approximate  formula  is 


2  J  +  2  £  ^  ] 

»nsf  U  yvx 


(2-33) 


The  signs  in  the  formulas  above  pertain  to  the  first  quadrant  in 
the  r\  -plane  (  t  >,  O  ,  6  >  O  ) ,  which  maps  into  a  rectangle  in  the  first  quad¬ 
rant  of  the  ^  -plane,  with  sides  located  on  the  lines 


*  K(4)  -  f  •  ■=  • - 

J0  VT - 1*  V>  -  S'*  t‘ 


~tf  ~j - i  f 

2  ZJ0  V;  -t*  Vi  ~  Ci  -S*)  t5 


(2-34) 


The  remaining  three  quadrants  in  the  rt -plane  map  into  the  remain- 
ing  three  quadrants  in  the  ^  -plane,  as  described  in  Reference  (13),  p.  377;  the 
cuts  from  t  S  to  the  unit  circle  along  the  real  axis  become  the  left  and 
right  sides  of  the  rectangle  in  the  ^  -plane: 


13.  Erdelyi,  A.,  et  al.  Higher  Transcendental  Functions,  Volume  2,  p.  377, 
McGraw-Hill  Book  Company,  New  York  (19S3). 
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Figure  13.  The  and  <  Planes 


Finally,  the  ^  -  plane  is  re-normalized,  so  as  to  lie  between  -1  and  +1  on 
both  axes : 

**«  '  KC*>  *  ‘  i  K\*t  c2-35) 


A  grid  is  now  to  be  set  up  in  the  £  plane,  and  mapped  back  to 
the  2  -  plane.  This  process  is  facilitated  by  first  rearranging  the 
quadrants  in  the  %  -  plane,  by  using  the  periodicity  of  the  elliptic  functions 

A  a/ 

as  follows:  in  the  first  and  second  quadrants,  let  %  -  £  ,  and  in  the  third 
and  fourth  let  $  =  -  (  f  *  2  K  (it)  )  and  use  the  relations  (see  for 
example,  Eqs.  122.00  and  122.04  of  Reference  11:) 
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5n  (  %)  •  -  sn(  ^  *  2  K )  =  S/i[-f^+aK)J  =  Sn  +  K)J  (2-36) 

A 

The  %  plane  then  has  the  form: 


Figure  14.  The  £  Plane 

A 

Two  types  of  grid  can  be  selected  in  the  %  plane:  a  rectangular  one 
(if  ISHEAR=0)  or  a  sheared  one  (if  ISHEAR=1).  The  latter  is  the  default. 

For  the  rectangular  grid,  equally  spaced  points  are  assigned,  ac¬ 
cording  to 

■  **<£)  =  (K-OAirt  ,  K-  7,2,...,  KMX  (2-37) 

a  *  J  K'dfc)-  (Z.-7)  ,  L  *  t,2,...,LMX 

where 

-  4K(4)  p  _  g-  xVil) 

“  KMX-1  ’  ~  LMX-1  C2_38) 

Becuase  the  mapping  is  conformal,  the  images  of  these  grid  lines  will  intersect 
at  right  angles  when  mapped  back  to  the  physical  plane. 

The  rectangular  grid  has  the  property  that  in  general  the  trailing 
edge  is  not  connected,  by  a  grid  line,  to  the  point  at  downstream  infinity. 
However,  such  a  connection  is  a  desirable  feature  in  certain  flowfield  codes 
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(see  Ref.  14,  for  example).  To  allow  this  feature,  the  ISHEAR=1  option 
establishes  a  sheared  grid,  consisting  of  the  same  lines  as  above,  but 

replacing  the  &  Cf)  lines  by  a  set  of  parabolas  which  intersect  the  blade 
surface  at  90  degrees,  and  are  displaced  from  a  base  parabola  that  connects 
the  trailing  edge  to  the  image  of  downstream  infinity: 


This  grid  is  given  by 


I)4i«  ,  t 


f  (TB)  t  (K-1)  A  f 


Const 


K  -  7,  Z,  ..  XMX 

L  -  7,2,...  ,  iMX 


(2-39) 


where 


- T  ;  -  &  |<r  [tfj  . ,]} 

4  [dK(4)  +  cr£j]  *  L  J" 


'(2-40) 


Some  points  on  these  parabolas  will  lie  outside  the  range 


-3K(*)  *  $  +K<&) 

When  this  occurs,  equivalent  points  are  found  by  adding  or  subtracting  4K(i), 
the  period  of  the  elliptic  sine.  In  addition,  the  base  parabola  is  always 
joined  to  the  lower-left  corner  of  Fig.  15,  by  subtracting  4K(4)from  the 
real  part  of  g  ,  if  the  latter  is  greater  than  -/<(■£). 

A 

This  completes  the  definition  of  the  grid  in  the  £  -plane.  Each 
of  these  grid  points  must  now  be  mapped  back  to  the  physical  plane.  The  first 


14.  Nenni,  J.P.  and  Rae,  W.J.,  "Experience  with  the  Development  of  an  Euler 
Code  for  Rotor  Rows",  ASME  Paper  83-GT-36  (March  1983). 
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transformation  is: 


(2-41) 


q  =  *sn(£  ;  i)  =  .TsnC?  ;  *) 

The  elliptic  sine  of  a  complex  argument  is  expressed  as  (Ref.  11,  Eq.  125.01) 

Sn  (u  +  4  tr,  •£) 

Sn  (u.,  £  )  <±n(v ,  &')  +  L  Cn(u.,-k)  aLn  (u,  -k)  sn  (  V,  Sk ')  CnCv,-k') 

1  -  S^Cv.-A')  oC^CU.-k)  (2-42) 

The  functions  in  this  expression  are  evaluated  by  the  Arithmetic-Geometric 
Mean  method  (see  Ref.  15,  p.  571): 

Set 

cl0  =  r  ,  ba  =  -k'  ,  ca  -  -A 

and  then  calculate 

a'ri  *  *2*  C^rt  -  i  +  ^  n  -  i  )  »  *  V<X0  _f  ' 

Cn  *  (2-43) 

until  Cn  =  0  to  a  prescribed  tolerance  (io~7  was  used  in  the  present  case.) 
Then  form 

n  -  2  •  An  •  u. 


and  calculate  >  <PM 

-  2  >  •“)  9o 

from 

<9  =  ±  « 

2 

^fn  +  ay'jC^*yr>  ( 

(2-44) 

Then  the  desired  results  are  given  by 


15. 


Abramowitz,  M.  and  Stegun,  I. A.,  Handbook  of  Mathematical  Functions, 
National  Bureau  of  Standards,  Applied  Mathematics  Series  55  (1964). 


sn(u,&)  -  (fo  t  cntu.-k)  =  C&a  cf0  } 


oLn  (ll,  -k  )  = 


<**  <fo 

(<9,  -Q.) 


(2-45) 


These  values  of  r\  are  then  mapped  back  to  the  i,  -  plane 


/S'?  -  a-*' 


1  "  * 


and  thence  to  the  uJ-  plane  by 


to  *  Z,  2  ( i  *  i-  B  ■  )  £  * 

~  A  m  n 


Lastly,  the  value  of  2  must  be  found  by  inverting: 


(2-46) 


(2-47) 


U  =  C 


uj-a. 
co  -  b 


-4^  TT 


-44-m  rr 


H  +  iG 

z  -z* 

H  +  i& 


3  C$(2>]  K  (2-48) 


This  is  written  as 

/r( z )  =  ^cz)  -  jo/ 

and  is  solved  by  a  Newton- Raphson  procedure 


where 


oo  *»; 


r 


(n;  r(2tn,) 

Z  "  F'Cz'”’) 


1t(iT  —  Hw) 

- - - 

7j-  H  +  l  Gr 

H  +  LG  .  a  _  2 -2* 

7T  - 

H  +  L  Gr 


(2-49) 


(2-50) 


(2-51) 


The  final  result  of  this  process  is  shown  in  Fig.  15,  for  a 
(examples  of  sheared  grids  can  be  found  in  Refs.  4  and  15). 
grid  line  in  Fig.  15  which  goes  to  downstream  infinity  does 
at  the  trailing  edge. 


rectangular  grid 
Note  that  the 
not  originate 


Figure  15.  Coordinate  Mapping 


:)i 
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Section  3 

METRIC  EVALUATIONS 


The  coordinate  transformation  enters  the  flowfield  solution  algorithm 
only  in  the  metric  derivatives.  These  can  be  evaluated  by  differencing  the 
coordinates  themselves,  or  in  the  case  of  a  conformal  transformation,  by  evalu¬ 
ating  the  analytic  expressions  for  them.  These  analytic  expressions  are  derived 
in  Ref.  1,  and  the  code  listed  in  that  report  contains  the  Fortran  statements 
required  to  evaluate  these  expressions.  However,  as  pointed  out  in  Ref.  14, 
the  truncation  error  resulting  from  the  use  of  analytic  metrics  in  the  finite- 
difference  flowfield  code  is  large  enough  to  cause  major  instabilities  in  the 
solution  algorithm.  It  is  highly  preferable  to  use  metrics  that  are  found  by 
differencing  the  coordinates  in  the  same  manner  as  the  flowfield  variables  are 
differenced.  A  program  to  achieve  this,  for  the  grid  conventions  used  in 
Ref.  14,  is  given  in  Appendix  D. 
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Section  4 
COMPUTER  PROGRAM 


A  listing  of  the  computer  program  is  given  in  Appendix  B,  and  a 
dictionary  of  variables  is  given  in  Appendix  C.  This  section  contains  a 
general  description  of  the  program,  plus  some  specific  details. 

In  order  to  handle  complex  arithmetric,  all  variables  beginning 
with  the  letter  Z  are  declared  to  be  complex  by  an  implicit  type  specification 
at  the  beginning  of  the  program. 

The  input  is  generally  described  by  comment  cards  in  the  deck.  Cer¬ 
tain  blade-shape  parameters  must  be  read  in:  EX,  G,  H,  ZLE,  ZTE,  ZN,  ZT.  Also, 
if  IGOT=l,  ZC  must  be  read  in.  The  blade  shape  itself  is  defined  by  pairs  of 
coordinates  in  the  S  ,  n.  plane  -  KJS  on  the  suction  side  and  KJP  on  the  pressure 
side.  In  the  version  shown  here,  these  were  read  in  as  a  table  in  subroutine 
SHAPE. 

The  blade  surface  coordinates  are  numbered  from  1  to  KJMX;  point  number 
1  is  the  trailing  edge,  2  through  KJIM  are  on  the  pressure  side  from  trailing 
edge  to  leading  edge,  point  KJLE  is  the  leading-edge  point  (ZLE),  KJLP  through 
KJMXM  are  on  the  s  jtion  side  from  leading  edge  to  trailing  edge,  and  point 
KJMX  repeats  the  trailing  edge. 

The  complex  sine  functions  are  calculated  next;  then  the  iterations 
to  determine  the  parameter  C  are  done.  The  initial  guess  provided  for  C  is 
He  =  f 4-t  .  it  may  happen  in  som.»  cases  that  a  better  guess  is  required: 
in  particular,  it  is  necessary  that  the  value  of C  must  lie  outside  the  blade- 
surface  curve  in  theJ2.  -plane.  If  it  does  not,  then  the  interior  of  this 
curve  in  thefl  -plane  is  mapped  to  the  exterior  of  the  blade-surface  image  in 
the  u>  -plane.  This  fact  can  be  seen  from  the  discussion  by  Kober  (Ref.  16) 
of  the  bilinear  transformation  applied  to  circles. 

16.  Kober,  H.,  Dictionary  of  Conformal  Transformations,  Dover  Publications, 
New  York  (1957 T- 
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After  the  parameter  c  has  been  found,  the  transformation  to  the  £ 
plane  is  carried  out,  using  the  fast  Fourier  transfrorm  procedure  (see 
Appendix  A) .  The  actual  calculations  of  the  Fourier  coefficients  are  done 
in  subroutine  FFT2,  a  proprietary  program  of  International  Mathematical  and 

Statistical  Libraries,  Inc.  (IMSL).  This  routine  computes  the  fast  Fourier 
transform  of  a  complex  vector  of  length  equal  to  a  power  of  two  (here  26) . 

The  coefficients  of  the  input  vector  are  given  in  normal  order  by  the 
array  named  as  the  first  argument  of  the  call;  the  coefficients  of  the 
output  vector  are  overstored  in  this  array,  in  reverse  binary  order. 

The  subroutine  SHUFL  is  then  used  to  restore  this  output  to  the  normal  order. 
The  coefficients  in  the  series  expression  for  ^  are  determined  iteratively 
in  a  relaxation  process  that  is  terminated  when  the  maximum  change  in  &  falls 
below  the  tolerance  ANGERR,  or  when  IMX  iterations  are  done. 


In  doing  these  iterations,  it  is  necessary  to  know  values  of 
in*  r  at  given  values  of  <9" ;  these  are  found  by  a  spline  fit  in  subroutine 

CISPLN,  which  is  a  straightforward  implementation  of  the  formulas  given  by 
17 

Ahlberg,  et  al. 


In  certain  cases  (typically  when  RMAX/RMIN  is  large)  the  calculated 
variation  of  O  with  may  be  non-monotonic;  if  this  occurs,  the  calculation 
should  be  repeated,  with  a  smaller  relaxation  factor.  The  progress  of  the 
iterations  is  printed,  showing  at  each  iteration  the  largest  change  in 
9  and  the  number  of  reversals  (i.e.,  the  number  of  occurrences  of  non-monotonic 
variation) . 


Next,  the  parameters  and  are  found,  starting  with  "best-guess" 
values  calculated  in  the  sectors  described  in  Section  2.  Once  these  are 
found,  the  mapping  to  the  fj,  -plane  follows.  The  calculations  that  link  the 
XL  -plane  and  the  £  -plane  are  done  in  subroutine  OMETA,  which  sums  the 
Theodorsen-Garrick  series,  using  complex  arithmetric.  This  subroutine  has 
been  modified  slightly  from  that  appearing  in  Ref.  1,  as  follows:  the  previous 
code  evaluated  sums  of  the  form 


17.  Ahlberg,  J.H.,  Nilson,  E.N.,  and  Walsh,  J.L.,  The  Theop'  of  Splines 
and  Their  Applications  Academic  Press,  New  York  (1967) . 
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2  SUM  =  £  zee  (JP)  z; 

Jp =  1 


JP  -  J 


(4-1) 


by  a  sequence  of  multiplications  and  additions: 

2  SUM  =  Zee  U>S)  ^  +  Zee  C6>4)  %U3+  z  cc  (1  )  (. 

The  current  version  uses  (Ref.  18,  p.  28) 

2.  SUM  *  <  *j[?CC  ^5KtZCC  (<^4)  J  £  +  Z  CC(65)j»  £  +  ... 


(4-3) 


Finally,  the  blade-surface  image  is  mapped  into  the  %  plane,  using 
the  elliptic-function  formulas  of  Nielsen  and  Perkins10  and  of  Luke,1^  as 
outlined  in  Section  2.  This  completes  the  mapping  of  the  blade  surface. 


It  is  now  possible  to  set  up  a  grid  in  the  f  -  plane,  and  map  it 
back  into  the  physical  plane.  This  involves  straightforward  evaluations  of 
the  transformation  functions.  The  only  new  complications  are  the  need  to 
evaulate  the  Jacobian  elliptic  sine  (done  in  subroutine  JCELFN)  and  to  provide 
a  good  initial  guess  for  the  Newton-Raphson  iteration  used  in  finding  2(i1)  . 
This  guess  is  provided  by  starting  each  series  of  calculations  on  the  image 
of  the  blade  surface  in  the  £  plane,  and  interpolating  to  find  the  corresponding 
point  in  the  S,  n.  plane.  The  interpolation  is  done  as  follows:  the  blade- 

A 

surface  image  values  of  %  are  stored  in  the  array  ZA1(I),  where  I  ranges  from 
1  to  KJMX.  Each  of  these  maps  into  a  point  in  the  S  ,  n  plane,  which  is  stored 

A 

in  the  array  ZA2  (KJ),  where  KJ  also  ranges  from  1  to  KJMX.  Thus  the  £  plane 
has  the  appearance: 

I  a  KSDCE-U 


I  =2 

J  I  =  KJMW 

rDXIR 

I.-lf  KJMX 

h-*i  f 

x  *  Keoee 

DXIM 

.L  =  i- MX 


Hartree,  D.R.,  Numerical  Analysis,  2nd  Edition,  Oxford,  Clarendon 
Press  (1958) . 


Points  in  the  uniform  grid  are  denoted  by  K,  L,  where 


(U(£)  =  -  3K  <&)  +  (K-  1)  DXIR  1,  KMX 

A  K  (4t  ) 

JU(£;  =  — —  -  (L-1)  DX1M  L=  1,  LMX. 


(4-4) 


The  value  of  I  for  which  the  real  part  of  ZA1(I)  is  largest  is  denoted  as 
KEDGE.  Then,  for  a  chosen  value  of  K  (with  L  =  1)  the  array  ZA1(I)  is 
searched  (first  for  I  =  KEDGE  to  KJMX,  then  for  I  =  2  to  KEDGE- 1)  until  a 
value,  called  I  ,  is  found  such  that  <8e[2-0KI)  <  &  (%  )]  .  if  none  is  found, 

*  -k 

it  is  concluded  that  I  =  KEDGE.  Linear  interpolation  is  then  used,  between 

■k  k 

the  points  I  and  I  -  1,  to  provide  the  initial  guess.  Fifty  iterations 
are  allowed  for  this  step,  at  each  value  of  K  and  L.  If  no  solution  is  found, 
the  value  (0,0)  is  printed. 


The  calculation  of  the  images  of  the  grid  points  in  the  various 
planes  is  bypassed  for  the  points  at  upstream  and  downstream  infinity,  and  at 
the  trailing  edge.  (The  trailing-edge  point  will  be  a  grid  point  if  ISHEAR=1). 
The  3L  -plane  locations  of  upstream  and  downstream  infinity  are  arbitrarily 
assigned  to  finite  locations  given  by  linear  extrapolation  from  the  two  adjacent 
L  -values. 


The  point  at  upstream  infinity  will  be  a  grid  point  only  if  KMX  is 
odd;  in  this  case  IOE=l,  and  the  image  calculations  are  bypassed. 


Adjustments  to  the  grid-point  image  locations  in  the  2  -plane  are 
sometimes  required  for  points  on  the  periodic  ouundary  (L=LMX)  for  values  of  K 
near  1,  KMX/2,  and  KMX.  At  these  points,  the  Newton-Raphson  iterations  used 
in  going  from  the  fl  -  plane  to  the  £  -plane  sometimes  cannot  distinguish  between 
points  that  are  separated  by  H  +  l  ( r  .  The  problem  can  be  seen  best  in  the 
U)  -plane.  (The  sketch  below  is  for  an  even  value  of  KMX;  the  same  picture 
applies  for  an  odd  value): 


In  Ref.  1,  a  quadratic  interpolation  was  used.  This  introduced  many  coiqplica- 
tions,  in  order  to  avoid  conditions  where  the  interpolation  base  points 
straddle  a  sharp  leading  or  trailing  edge,  and  where  I*  is  near  KEDGE  or 
KEDGE- 1.  The  logic  in  Ref.  1  was  not  adequate  to  handle  all  such  cases. 
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Note  that  the  same  point  in  the  uj  -plane  (and  thus  also  in  the H -plane)  can 
map  into  either  of  two  points  in  the  £  -plane,  which  differ  by  H  +  i  &  .  The 
selection  is  guided  toward  the  correct  value  by  starting  on  the  blade  (L=l) 
and  working  toward  the  periodic  boundary  (L=LMX) ,  but  it  can  happen  that  the 
wrong  branch  is  chosen  during  the  iterations.  To  avoid  this,  the  imaginary 
parts  of  2  and  IN  are  compared,  for  K  values  near  the  leading  edge,  and  the 
imaginary  parts  of  £  and  £T  are  compared  near  the  trailing  edge,  and  the 
quantity  t-55is  added  or  subtracted  (depending  on  the  value  of  ^  )  where 
necessary. 

Finally,  the  real  and  imaginary  parts  of  Z  are  punched  on  cards 
(if  PNCHZA=TRUE) .  These  values  can  be  used,  in  a  separate  program,  to  calculate 
the  metrics  of  the  transformation  (see  Appendix  D) . 
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4. 


Section  5 

CONCLUDING  REMARKS 


The  program  described  above  has  been  applied  to  very  few  cases  - 
the  thick  blades  used  for  illustrative  purposes  in  this  report,  and  the  cases 
shown  in  Refs.  4  and  15.  Because  of  this  limited  experience,  the  range  of 
applicability  of  the  program  is  largely  unknown.  On  the  basis  of  current 
experience,  it  appears  that  the  Theodorsen-Garrick  step  may  not  converge  for 
gap/chord  ratios  less  than  around  0.8.  The  set  of  numerical  tolerances, 
maximum  iteration  counts,  and  relaxation  factors  used  may  need  to  be  adjusted 
for  certain  shapes.  In  addition,  new  approximations  may  be  required  as  initial 
guesses  in  various  places,  in  order  to  handle  such  features  as  leading  edges 
with  very  small  radii  of  curvature. 


I 


Appendix  A 

DETAILS  OF  THE  FAST  FOURIER  TRANSFORM  PROCEDURES 


Equations  2-20  and  2-21  contain  ZN+2  constants,  which  are  evaluated 
as  follows:  first,  they  are  satisfied  at  a  discrete  number  of  points,  denoted 

by  : 


K  - 


znCK-  i) 
~2N 


K 


ZN 


(A-1) 


where N  was  chosen  to  be  64,  in  the  present  case.  Thus: 

K 

rK  =  R0  +  Z  (A;  1  4>k  -  6;  J  $  )  +(-i)  RH  (a-2) 

r1 

I  i  ** +  ^  i**'*  +(-nK  bn  (a.3) 


Each  right-hand  side  now  contains  2  N  coefficients.  The  correspondence  between 
either  of  these  and  the  Fast  Fourier  Transform  (FFT)  as  presented  in  Ref.  6 
is  given  by  the  next  two  equations:  consider  the  expression 


2N- 1  nj 

Y  <1>  -  L  C  in)  W2n 

n  «  o 


i 


*  0,1  y  2, 


2  /V  -  7  ;  W 


IN 


Ztt 

2.N 

(A- 4) 


where  values  V  (• )  are  real,  and  the  2.N  values  of  C(0  are  in  general  complex, 
but  must  satisfy  the  following  redundancy  condition,  in  order  that  the  y  <*• 
values  be  real: 

Cin )  =  c(ZN-n)  ,  n  =  7,2,...,  A/-;  (A-5) 


C  <o)  and  C-Of)  pure  real 

where  the  tilde  denotes  the  complex  conjugate.  When  these  conditions  are 

met,  Eq.  A-4  can  be  written  as 

I  N-i  r  tiirc  nlTC 

Yfi;  *  CRio)  +(->)  Cs<N)  +  2£  «jC*(n)C<w  — - C^n) 

i«  <7,  7,...,  2A/-J  (A-6) 

This  form  can  now  be  used,  in  conjunction  with  Eq.  2-20  or  2-21,  to  facilitate 
application  of  the  FFT  to  the  complex  form  given  in  Eq.  A-4. 


A-1 


In  the  case  of  Eq.  A-2,  values  of  the  coefficients  R0  through 
and  6,  through  6w_t  are  found,  from  given  values  of  r  .  This  is  procedure 
4  of  Ref.  6,  which  takes  the  following  steps: 


1.  Set  X/*)  -  V(2&)  =  (Jt~'rKk 

X,(4)  =  y  (2-feft)  -  <^r>2*  +  f 


Je  *  o,7,  ...,  A/  -  7 


2.  Set  =  Xl(f)+iXz  (j)  ,  ^  =  <9,  7 ,  .  .  N  -  7  (A-8) 

3.  Calculate  the  W  -point  Discrete  Fourier  Transform  of 

•  1  N~1  ~n4 

fRcn)  -  +  i/?)z(n)  -  —  JT  X(j)  WN  ;  «  «  O,  /V  -  J 

.  L  ^ 

WN  = 

4.  By  periodicity,  set  ^  f A/ >  *  /f?  (tf) 

5.  Apply  Eq.  34  of  Ref.  6,  in  order  to  extract  ^?f(«)and  ^)  (*»)  from 

.  r  1  *  *■  *  0,1,  (A-10) 

flz(n)  =  j  -j/?  (A /-  r)  -  (a;| 

Note  that  these  expressions  use  ^C«?for  n  *  0,1,..,  N  to  give  >#f(A)  and  2(a) 
for  n  •  <?,  1 ,  ... }  N/z 

6.  These  values  of  /$,(»)  and  A^n)  then  give  C(.n)  for  the  same  range: 

C(A>  -  Y  ,  A  *  0,1,...,  N/a  (A- 11) 

N 

For  the  range  of  A  from  +  7  to  A/—  7  ,  use  Eq.  36  of  Ref.  6,  with  A  replaced 

* 

by  M-  A  (and  noting  that  W  «  -  1  ): 


C  (n)  -  C(2A/-a) 


(A- 12) 


7  a  -  4"  +  1 »  7~  +  2/---^-; 


A/ 

This  equation,  applied  for  rt  *  J-+1  , -g +2  uses  |P,  and  fi2  with  index 

E-1 ,  1  to  get  C(n. )  for  n*  ^*1,  ,  N~1  .  The  process  is  completed 

by  setting  r  _  1  (A-1J) 


7.  The  Rj  and  coefficients  are  retrieved  from: 

fla  *  [c  «»>]  ,  «  &.  [c(*>j  ^  ^ 

Rj  =  z(U[c^)]  ,  B-  =  2  [c  Cj>]  ,  Jt2,...,N-l 

At  this  point,  Ba  and  are  undetermined.  Following  Ives,  B#  is  set  equal 
to  zero,  and  B0  is  chosen  so  as  to  place  the  trailing  edge  at  <p*0  (this  latter 
selection  of  B0  is  actually  carried  out  in  a  subsequent  step,  noted  below). 


In  the  case  of  Eq.  A-3,  the  /? ' s  and  B  's  are  considered  known,  and 
are  used  to  evaluate  .  The  coefficient  &a  can  be  found  from 


Ort  *  +  2)  ( &N  -  ° ) 

i*1 

Actually,  it  is  simpler  to  evaluate  the  right-hand  side  of 


(A- 15) 


^-k'Ba  =  ^  i  +  1  4+  )  >  4  •  N-1 

ja1  (A- 16) 

and  then  find  B„  from 


0.  -  0.  -  bus). 


(A- 17) 


The  actual  ^valuation  of  the  right-hand  side  takes  the  following  steps  (Pro¬ 
cedure  5  of  Ref.  6):  by  comparison  of  Eqs.  A-3  and  A-6: 

1.  Set  C  Co)  -  O ,  CCM)  =  O 

i  r_  •  -  i  .  .  ..  (A- 18) 


n>  s  2  [Bn  -  1  qn  ]  >  n  •  1,2,...,  N 


Note  that  the  C  '  s  determined  here  are  different  from  those  used  in  Procedure  4; 

the  flj  and  S.  values  are  the  same,  but  their  relation  to  C.  is  different. 

1  <  i 


2.  Values  of  are  then  used  to  find  and  fi2Cn)  :  Equations  40  and  41 

nf  Ref.  6  are  rewritten,  using 


C(/t)  *  c(2  A»-/t) 


Replace  n  by  N  +  n  : 


C  ( N  +n)  *  c(ZN-N-n)  =  c(N-r t) 


Thus  Eqs.  40  and  41  of  Ref.  6  are 

=  C(n)  +  c  ( N -n) 

=  [c(*o  -  C  (A/  -  *  )J  W2N 
These  are  all  the  values  needed  for  and  . 


n  l  N/2 


(A-19) 


3.  Find  /R(n)  ,  n.~ori from  Eqs.  42  and  43  of  Ref.  6: 


fl(N-n)  =  £,(.*)  +  i  fi£(n) 


n  *  0,  /,-•'>  A//2 


This  gives,  on  the  left-hand  side,  all  values  from  rt=o  to  «  »A/ . 


4.  Calculate 


(A- 20) 


X(j>  *  Z  ^  Y  -  O }  A/  -  T 


(A- 21) 


5.  Finally: 

~  Ba  •  (L  [XC«] 

►  -ft  »  0 1  i f  ... ,  /V -r 

-  *.*♦,  '  s.  -  [xrt>]  J  <*-22> 

The  first  of  these  equations,  with  \*o,  is  used  to  find  Ba  . 


In  order  to  apply  these  formulas,  it  is  necessary  to  have  a  relation 


A- 4 


between  f  and  &  : 


rK  *  r(°K^  >  X  '  1,2,  ■■■>  ZH  (A- 23) 

which  is  found  from  a  spline  fit  to  the  blade-surface  image  in  the  uJ -  plane. 

The  discrete  Fourier  transform  and  its  inverse  are  given  by 
;  /V  -  J  - 

X^>  WN  ,  rf  0,1 ,2,...,  N-1  (A- 24) 

X^>  -  £  WN*  ,  i  •  O,  1, 2,...,  N-  1  (A- 25) 

n  s-o  ® 

The  IMSL  routine  FFT2  evaluates  the  second  of  these,  i.e.,  it  returns 
X(^)  ,  given  the  values  /Pin.’)  .  To  evaluate  the  first  of  these,  FFT2  is 
used  with  input  X(^)  /  H  ,  and  with  output  interpreted  to  be  (R(n)  ; 
this  is  Procedure  1  of  Reference  6: 


-v  ~ 

NSHn)  -  £  X(f)  WN 


(A- 26) 


In  the  FORTRAN  version  of  these  and  other  procedures,  it  is  conve¬ 
nient  to  use  indices  that  begin  at  one,  rather  than  zero,  by  setting  =  ^ 
(FORTRAN  symbol  JP).  The  corresponding  table,  for  example  of  the  coefficient 
B. ,  is 


A- 5 


j 

JP 

B. 

1 

B(JP) 

0 

1 

B 

0 

BCD 

1 

2 

B1 

B(2) 

2 

3 

B2 

B(3) 

N-l 

N 

bn-i 

B(N) 

N 

N+l 

bn 

B(N+1) 

In  addition,  care  must  be  taken  with  the  argument  N~n;  for  example,  Eqs. 
are  written  as 


£  ft  1  ( NP  ) 

ZftZ(NP) 


2  CC  (NP)  +  ZCC(N  +  2~Np) 
[ 2CC(NP )  -  ZCC  (N+Z-  NP 


NP-  1,2,..., 


(A- 


It  can  be  verified  that  the  quantity  N  +  2  -  NP  in  the  arguments  above  pre 
serves  the  correct  ordering  -  for  example  when  N  =  64,  and  n  =  0,  C,(N~n ) 
This  would  be  stored  in  ZCC (64  +  2  -  1)  =  ZCC (65). 


A- 6 


1 


27) 

C^4) 


Appendix  B 

COMPUTER  PROGRAM  LISTING 


LEVEL  21.7  (  DEC  72  ) 


.ER  OPTIONS 


OS/360  FORTRAN  H 

NAME-  MAIN,OPT-02,LINECNT-60,SIZE-0000K, 

SOURCE .EBCDIC, NOL 1ST, NODECK, LOAD, MAP .NOEDIT, ID, XREF 


ISN  0002 
ISN  0003 
ISN  0004 
ISN  0005 
ISN  0006 
ISN  0007 
ISN  0003 
ISN  0009 
ISN  0010 
ISN  0011 


ISM  0012 
ISN  0013 
ISN  0014 
ISN  0015 

ISN  0016 


ISN  0017 
ISN  0013 
ISN  0019 
ISN  0020 
ISN  0021 
ISN  0022 
ISN  0023 
ISN  0024 
ISN  0025 
ISN  0026 
ISN  0027 
ISN  0029 
ISN  0023 
ISN  0030 
ISN  0031 
ISN  0032 


PROGRAM  IVLTMP  -  THE  IVES  -  LIUTERMOZA  CONFORMAL  TRANSFORMATION  FOR 
TURBOMACHINERY  CASCADES  (AIAA  JOURNAL, VOL.  5,1977,  PP  647  -  652) 
DOCUMENTATION  IS  GIVEN  IN: 

W.  J.  RAE,  A  COMPUTER  PROGRAM  FOR  THE  IVES  TRANSFORMATION 
IN  TURBOMACHINERY  CASCADES,  CALSPAN  CORPORATION  REPORT  6275-A-3, 
NOVEMBER  1980 

W.J.  RAE,  MODIFICATIONS  OF  THE  IVES  -  LIUTERMOZA  CONFORMAL- 
MAPPING  PROCEDURE  FOR  TURBOMACHINERY  CASCADES,  ASME  PAPER 
83-GT-U6,  MARCH  1983 

W.J.  RAE,  REVISED  COMPUTER  PROGRAM  FOR  EVALUATING  THE  IVES 
TRANSFORMATION  IN  TURBOMACHINERY  CASCADES,  CALSPAN  CORPORATION 
REPORT  7177-A-l,  JULY  1983 

IMPLICIT  REAL*8< A-H ,0-Y  > , COMPLEX* 16( Z ) 

LOGICAL  PNCHZA 

C0MM0N/TGINTG/N.NP1 , NP2 . N2 , NB2 , NB2P 1 , IP, IPMX, ITP . IWK, IMX.KJMX 
C0MM0N/TGCMPX/Z1 .ZW2N.ZI , ZNN , ZA , ZCC , ZA1 ,ZA2 
COMMON /TGDBL E /P BN , OM ,QMM, ANGERR , A ,B,E,F,THT,PHI ,X,Y 
DIMENSION  ZS(80),ZP(80), ZOMS ( 80 ) , ZOMP  (  80 ) 

DIMENSION  RDS(80),RDP(80),THS(80).THP(80) 

DIMENSION  RDSX( 80  > , RDPX<  80  > ,THSX( 80  > ,THPX<  80  > 

DIMENSION  X< 160),Y( 160),E( 1 50 > ,  F ( 1 50 > , THT< 1 60 > 

DIMENSION  PH  I  (  130),A<65)  , B< 65 > ,ZCC< 65 > , 

*  ZA< I65),ZA1< 165>,ZAZ(  165), 

*  IWK<  7 ) 

DIMENSION  ITP< 100  ) , ID( 36 ) 

DIMENSION  XX(50,20),YY(50,20) 

DO  10  I  -  1,7 

10  IWK< I )  -  0 

NAMELIST/ INPUTS/  ANGERR , EX , G , H , IGOT , ILE , ITE , KMX , 

*  IMX ,  LMX ,  OM,  PNCHZA,  RTOL ,  ZC,  ZLE ,  ZN,  ZT.  ZTE , 

*  IPMX, ISHEAR, INR.KJS.KJP.KNR 

—  SET  NAMELIST  DEFAULT  VALUES 

ANGERR-0 .0100 
IGOT-O 
ILE  -  1 
ITE-0 
KMX  -  40 
LMX  -  10 
KJS  -  20 
KJP  -  20 
IMX  -  400 
ISHEAR  •  1 
IPMX  -  1 
KNR  -  0 
INR  -  0 
OM-O. IDO 
PNCHZA-. FALSE. 

RTOL -3. 090 

NOTE:  EX.  G,  H,  ZC ,  ZLE,  ZN ,  ZT,  ZTE  HAVE  NO  DEFAULT  VALUES 
(ZC  IS  NOT  NEEDED  IF  IGOT-O) 


» 


ISM  0033 
ISN  0034 

ISN  0035 
ISN  0036 
ISN  0037 
ISN  0033 
ISN  0039 
ISN  0040 
ISN  0041 


ISN  0042 
ISN  0043 
ISN  0044 
ISN  0046 
ISN  0047 


ISN  0048 
ISN  0043 
ISN  0050 
ISN  0052 
ISN  0053 
ISN  0054 
ISN  0055 
ISN  0056 


C 

READ! 5 , 1 03 ) ( I D( I >,1-1,36) 

103  FORMAT! 18A4) 

C 

PI  -  4  .  ODO*DATAN(  1.000  ) 

TP  I  -  2 . ODO*P  I 

ZPI  -  DCMPLX! PI ,0.000 > 

Z1  »  DCMPLX! 1.000,0.000) 

ZERO  -  DCMPLXIO. 000, 0.000) 

ZMGA  -  DCMPLX!+1.0D0, 0.000) 

ZMGS  -  DCMPLX(-l.ODO.O.ODO) 

C  ***  READ  PNCHZA-T  IF  ( ZA<  L  ) , L-l , LMX  >  IS  TO  BE  PUNCHED  FOR  ALL  VALUES 
C  ***  OF  K,  OTHERWISE  READ  PNCHZA-F 

C  IGOT  -  1  IF  THE  VALUE  OF  ZC  IS  KNOWN.  OTHERWISE,  IGOT  -  0. 

C  KMX  AND  LMX  ARE  THE  GRID  SIZES  IN  THE  FINAL  TRANSFORMED  PLANE. 

C  KJS  AND  KUP  ARE  THE  NUMBERS  OF  POINTS  ON  THE  SUCTION  AND  PRESSURE 
C  SIDES  AT  WHICH  PAIRS  OF  BLADE  COORDINATES  WILL  BE  INPUT. 

C  ILE  -  3  OR  1  FOR  A  SHARP  OR  ROUNDED  LEADING  EDGE,  RESPECTIVELV 

C  ITE  -  0  OR  1  FOR  A  SHARP  OR  ROUNDED  TRAILING  EDGE,  RESPECTIVELY. 

C  IMX  IS  THE  MAXIMUM  NUMBER  OF  ITERATIONS  ALLOWED  FOR  THE  PHI/THETA 
C  ITERATIONS 

C  INR  AND  KNR.NE.O  WILL  CAUSE  A  DIAGNOSTIC  OUTPUT  FOR  THE  FIRST  INR 
C  NEWTON-RAPHSON  ITERATIONS  ! TO  FIND  Z,  GIVEN  ZBOMK ) ,  AT  THE  STATION 
C  K  -  KNR.  THE  PROGRAM  WILL  THEN  STOP. 

C  IPMX.NE.l  CAN  BE  USED  TO  DISPLAY  THE  VALUES  OF  THETA  DURING  THE 
C  PHI/THETA  ITERATIONS.  AT  THE  ITERATION  NUMBERS  READ  INTO  THE 
C  ITP  ARRAY  BELOW. 

C  ISHEAR  -  0  GIVES  AN  ORTHOGONAL  GRID.  THE  LINES  K-l  AND  K-KMX,  WHICH 
C  START  AT  THE  TRAILING  EDGE,  DO  NOT  GO  TO  DOWNSTREAM  INFINITY. 

C  ISHEAR  -  1  SHEARS  THE  GRID  UNIFORMLY!  IN  THIS  CASE,  THE  K-l  AND  KMX 
C  LINES  DO  GO  TO  DOWNSTREAM  INFINITY. 

C  OM  IS  A  RELAXATION  FACTOR  USED  IN  THE  PHI/THETA  MAPPING.  USE  0.1, 

C  OR  A  SMALLER  VALUE  IF  THE  A  AND  B  ITERATIONS  FAIL  TO  CONVERGE. 

C  ANGERR  IS  THE  ANGULAR  TOLERANCE  (IN  RADIANS)  FOR  THE  PHI/THETA 
C  TRANSFORMATION.  A  REASONABLE  VALUE  IS  0.01 
C  RTOL  13  THE  TOLERANCE  FOR  THE  MAX/MIN  RADIUS  RATIO  IN  THE 
C  OMEGA  PLANE 
C 

C  —  FOR  A  SHARP  LEADING  EDGE  <  ILE-0  >  ZN  MUST  EQUAL  ZLE 
C  — -  FOR  A  SHARP  TRAILING  EDGE  (ITE-O)  ZT  MUST  EQUAL  ZTE 
C 

C  —  READ  NAMELIST  INPUT  DATA 
C 

REAO(S, INPUTS  ) 

ITP!  1  )  -  999 

IF! IPMX.NE . 1  )  REA0(5, 102 )( ITP! IP  ) , IP-1 , IPMX  ) 

102  FORMAT! 2014) 

IP  »  1 
C 
C 

KMXH-! KMX+ 1 >/2 
KMXL-KMXH 

IF ( MOD! KMX , 2  ) . NE . 0  )  KMXL-KMXH-1 
IOE  ■  MOO! KMX, 2) 

KAA  -  KMXH  -  2  -  IOE 
KBB  •  KMXH  ♦  3 
KMXM4  -  KMX  -  4 
KOUNT  -  0 


ISN  0057 
ISN  0058 
ISN  0059 
ISM  0060 
ISN  0061 
ISN  0062 
ISN  0063 
ISN  0064 
ISN  0065 
ISN  0063 
ISN  0067 
ISN  0063 


ISN  0069 
ISN  0070 
ISN  0071 
ISN  0072 
ISN  0073 

ISN  0074 
ISN  0075 
ISN  0076 
ISN  0077 
ISN  0078 
ISM  0079 
ISN  0080 
ISN  0081 
ISN  0082 
ISN  0083 
ISN  0084 
ISN  0085 

ISN  0086 
ISN  0087 
ISN  0088 
ISN  0089 
ISN  0090 
ISN  0091 
ISN  0092 
ISN  0093 
ISN  0094 
ISN  0095 

ISN  0096 
ISN  0097 
ISN  0096 


C 

C  INPUT  VARIABLES  EX,  G,  H,  ZLE ,  2N ,  ZT,  ZTE , 

C  THE  FORTRAN  STATEMENTS  IN  SUBROUTINE  SHAPE, 

C  AND  THE  VARIABLES  LISTED  IN  COMMON  BLOCK  GEOM  < IF  ONE  IS  BEING  USED) 
C  ARE  ALL  SPECIFIC  TO  THE  BLADE  SHAPE  BEING  USED. 

C 

C 

C  CALCULATION  OF  THE  BLADE  SHAPE 
C 

D  -  CDABS! ZTE-ZLE  ) 

CALL  SHAPE (D.H.G, EX, ZP.ZS.KJS.KJP ) 

SG  -  DSQRT <  H*H  +  G*G  ) 

SIZE  »  5 . 0D0*SG 

ZDA  -  DCMPLX!G/SG,-H/SG) 

ZGAMMA-DCONOG<  ZDA ) 

WR ITE ( 6 , 207  > 

207  FORMAT! 1H1 ) 

WRITE (6. 206  X I0< I >, 1-1 ,36) 

206  FORMAT! 30 X , 18A4  > 

WRI TE ! 6 , 240 )  D,H,G,EX,SG 

240  FORMAT! //I OX, ' BLADE-GEOMETRV  PARAMETERS  ARE!', 

*  //5X, 'ABS(ZTE-ZL£  )  - ' , 

*  F 1 0 . 5 ,  '  H  -  ‘ , F 1 0 . 5 ,  '  G  ■  ',F10.5,'  EX  - 

*  F  1 0 . 5 , '  SLANT  GAP  -’.FIO.S,//) 

WRITEI6, INPUTS) 

KJLE  -  KJP  +  2 
KJMX  -  KJLE  ♦  KJS  ♦  l 
WRITE ! 6 , 209 ) 

209  FORMAT!  8X , ' BLADE  COORDINATES:’, 

*  //9X , 'SUCTION  SIDE'  ,  /3X ,  ’ KJ ' ,7X, 'S* ,13X, ' N '  ,13X,  'X'  ,13X, ' V ' ,/) 
ZXY  -  ZGAMMA*ZLE 

WRITE ! 6 , 270  )  KJLE, ZLE, ZXV 

270  FORMAT !I5,1P4E14.5) 

DO  64  K  -  l.KJS 

KO  -  KJLE  +  K 
ZXY  -  ZGAMMA*ZS!K> 

WRITE! 6 , 270  >  KJ,ZS!K),ZXY 
64  CONTINUE 

ZXY  -  ZGAMMA'ZTE 

WR  ITE ! 6 , 270  )  KJMX, ZTE  ,ZXY 

WRITE  16,271 > 

271  FORMAT! //8X, 'PRESSURE  SIDE' ,/3X, 'KJ' ,7X, 'S' ,13X,  'N'  , 

*  1 3X  ,  '  X  '  ,  1 3X  ,  '  Y  '  ,  /  > 

ZXY  -  ZGAMMA*ZLE 

WR ITE ( 6 , 270  )  KJLE, ZLE, ZXY 
DO  61  K  ■  l.KJP 
KJ  -  KJLE  -  K 
ZXY  •  ZGAMMA*ZP!K) 

WR ITE ( 6 , 270 )  KJ,ZP<K),ZXY 
61  CONTINUE 
KJ  -  1 

ZXY  -  ZGAMMA*ZTE 
WR ITE ! 6 , 270  >  KJ, ZTE, ZXY 
C 

ZDN  -  DCMPLX! H ,G  ) 

ZZ  «  ! ZT-ZN  )/ZDN 
ZTN  -  ZPI*ZZ 


ISN  0099 
ISN  0100 
ISN  0101 
ISN  0102 
ISN  0103 
IS II  0104 
ISN  0105 
ISN  0106 
ISN  010.’ 
ISN  0108 
ISN  0109 
ISN  0110 
ISN  0111 
ISN  0112 
ISN  0113 
ISN  0114 
ISN  0115 

isn  one 

ISN  0117 
ISN  0118 
ISN  0119 
ISN  0120 
ISN  0121 
ISN  0122 
ISN  0123 
ISN  0124 
ISN  0 125 
ISN  0126 


ISM  0127 
ISN  0129 
ISN  0130 
ISN  0131 
ISN  0132 
ISN  0133 
ISN  0134 
ISN  0135 
ISN  0136 
ISM  0137 
ISN  0139 
ISN  0140 
ISN  0141 
ISN  0142 
ISN  0143 
ISN  0144 
ISN  0145 
ISN  0146 
ISN  0147 
ISN  0148 
ISN  0149 
ISN  0150 
ISN  0151 
ISN  0152 
ISN  0153 
ISN  0154 


CHI  -  PI*EX*(G*DREAL( ZT-ZN  >-H*D IMAG( ZT-ZN  >  )/(SG*SG) 

XA  *  PI*EX*(H*DREAL( ZT-ZN)  +G*DIMAG( ZT-ZN))/(SG*SG) 

R  *  DEXP( -CHI  ) 

ZPLUS  -  DCMPLX(R*DCOS(XA),-R*DSIN<XA> > 

R  -  l.ODO/R 

ZMINUS-  DCMPLX(R*DCOS(XA),R*DSIN(XA)> 

00  20  K  -  1 ,  KJ  S 
ZZT  -  ( ZS( K  )-ZT  )/ZDN 
ZETA1S  -  ZP I *ZZT 
ZETA2S  -  COS  I N( ZETA1 S  ) 

ZZN  »  (ZS(K)-ZN)/ZDN 
ZETA3S  -  ZP I *ZZN 
ZETA4S  -  CDSINf ZETA3S  ) 

ZFS  -  ZETA2S/ZETA4S 
RDS(K)  •  COABS(ZFS) 

THS(  K  )•  DATAN2( DIMAG(  ZFS  ) , DREAL ( ZFS  >  > 

20  CONTINUE 

DO  24  K  *  l.KJP 
ZZT  -  ( ZP  <  K  )  -ZT ) /ZDN 
ZETA1P  «  ZP I *ZZT 
ZETA2P  «  C0SIN(ZETA1P  ) 

ZZN  -  (  ZP ( K  > -ZN  )  /ZDN 
ZETA3P  -  ZP I *ZZN 
ZETA4P  -  CDS IN<  ZETA3P > 

ZFP  -  ZETA2P /ZETA4P 
RDP(K)  -  CDABS(ZFP) 

THP(K)-  DATAN21 D IMAG( ZFP  >, DREAL (ZFP  >  > 

24  CONTINUE 

C  NOW  ADD  THE  LEADING-  AND  TRAIL ING-EDGE  POINTS.  AND  STORE  THE  G(Z) 

C  ARRAY  AS  E( KJ  >*EXP( I*THT( KJ  >  > .WHERE  KJ=1,KJMX  AS  YOU  GO  FROM  TE  AROUND 
C  THE  PRESSURE  SIDE  TO  THE  LE  ( KJ»KJLE  )  AND  THEN  ALONG  THE  SUCTION  SIDE 
C  BACK  TO  THE  TE  AGAIN  <KJ»KJMX>. 

IF( ITE.EO. I  )  GO  TO  13 
E< 1  )”0 . 030 
THT( 1 J-O.ODO 
GO  TO  14 

13  ZETA2“CDSIN<ZPI*<  ZTE-ZT  )/ZDN  > 

ZETA4-CDSIN(ZPI*(ZTE-ZN>/ZDN> 

ZFP-ZETA2/ZETA4 

E< 1 )-CDA2S(ZFP  ) 

THT ( 1 )«DATAN2 (DIMAG( ZFP  >, DREAL (ZFP  >  > 

14  I F ( ILE.EQ.  1  )  GO  TO  15 
E(KJLE)  »  0.000- 
THT(KJLE)  -  O.ODO 

GO  TO  16 

15  ZETA2»CDSIN(ZPI*(ZLE-ZT)/ZDN> 

ZETA4*CDSIN(ZPI*(ZLE-ZN>/ZDN> 

ZFP-ZETA2/ZETA4 

E( KJLE  )  -  CDABS( ZFP  ) 

THT(KJLE)  ■  DATAN2( DIMAG( ZFP  ), DREAL ( ZFP  )  > 

16  DO  17  K  ■  l.KOP 
KJ  »  KJLE  -  K 

E ( KJ  )-ROP(K> 

THT(KJ)-THP(K) 

17  CONTINUE 

DO  18  K  ■  l.KJS 
KJ  »  KJLE  +  K 
E(KJ)-RDS(K> 


B-4 


I 


I SN  0155 
ISN  0155 
ISM  0157 


ISN  0158 
ISN  0160 


ISN  0162 
ISN  0163 
ISN  0164 
ISN  0165 
ISN  0167 
ISN  0168 
ISN  0169 


13  THT( KJ )«THS<  K ) 

E(KOMX)  «  E< 1 > 

THT <  KJMX  i  -  THT< 1 ) 

NOW  ADJUST  THE  BRANCHES  OF  G(Z>  SO  AS  TO  BE  CONTINUOUS  ACROSS 
THE  CUT  ( ALONG  THE  NEGATIVE  REAL  AXIS)  OF  THE  DATAN2  FUNCTION. 

IF<  H.LT.O.ODO  )  GO  TO  251 
I F <  ITE.EQ. 1 >  GO  TO  251 

FOR  COMPRESSORS  WITH  A  SHARP  T.  E.  ,  THE  CONVENTION  IS  TO  FORCE 
ARG<  THT  <  KJ "2  )  )  TO  BE  NEGATIVE: 

BR  =  0.000 
KA  =  3 
PO  -  THT <  2  ) 

IF(PO.LT.O.ODO)  GO  TO  252 
BR  =  -l.ODO 
THT <  2  )  -  THT <  2 )  -  TP  I 
GO  TO  252 


ISN  0170 
ISN  0171 
ISN  0172 
ISN  0174 
ISN  0175 
ISN  0175 
ISN  0177 
ISN  0178 
ISN  0180 
ISN  0182 
ISN  0184 
ISN  0185 
ISN  0185 
ISN  0187 
ISN  0188 
ISN  0189 
ISN  0190 

ISN  0191 
ISN  0192 
ISN  0193 
ISN  0194 
ISN  0195 
ISN  0196 
ISN  0197 
ISN  0198 
ISN  0199 
ISN  0200 
ISN  0201 
ISN  0202 
ISN  0203 
ISN  0204 


CONVENTION  FOR  ALL  TURBINES,  AND  FOR  COMPRESSORS 
WITH  A  ROUND  T.  E. ,  IS: 

251  BR=>0 .  ODO 
KA-3 

I F ( ITE . EQ . 1 >  KA-2 
PO«THT <  KA- 1 ) 

252  DO  321  KJ-KA.KOMX 
CHG-THTIXJ  >-P0 
P0»THT<KJ ) 

IF(DABS(CHG).LE.PI )  GO  TO  321 
I F ( CHG . GT . P I )  BR-BR-l.ODO 
IFfCHG.LT.-PI )  BR*BR* l . ODO 
321  THT!KO)»THT!KJ)+BR*TPI 
DO  323  K« 1 , KJP 


THP(K)  - 
323  CONTINUE 
'  DO  322  K 
THS(K)  - 
322  CONTINUE 


-  THT<  KJLE-K ) 


-  1  ,  KJ  S 
THT( KOLE+K  > 


DO  25  K  -  l.KJS 
ARS  »  EX*THS<  K  > 

RS  «  RDS( K  )**EX 
RDSXIK)  ■  RS 
THSX(K)  »  ARS 

25  CONTINUE 

00  26  K  *  l.KOP 
ARP  -  EX*THP<  K  ) 

RP  -  RDP ( K  >  **EX 
RDPX(K)  «  RP 
THPX(K)  -  ARP 

26  CONTINUE 
WRITE! 6,241 > 

241  FORMAT! //10X, ' BLADE -SURFACE  IMAGES  IN  THE  G  -  PLANE  IRATIO  OF' 

*  •  SINES)  AND  CAP  OMEGA  PLANE ! G** 1 /KAPPA >  .AND', 

*  //  9X , '  RADII  AND  ANGLES  USED  IN  SELECTING  THE  PROPER', 


1 


ISN  0250 
ISN  0251 
ISN  0252 
ISN  0253 


'  BRANCHES  OF  THE  RATIO  OF  SINE  FUNCTIONS  ARE : ' . 

//3X, 'KJ ' . 15X, 'G ' ,25X, 'CAP  OMEGA' . 13X, 'R' , 1IX, 'THETA' , 


*  9X, ' R**EX ‘ ,7X, ‘ THETA*EX ' , ) 

ISN 

0203 

XS-E<  1 )*DCOS! THT<  1)) 

ISN 

0206 

VS-E<  1 ) *DS IN ( THT(  1>> 

ISN 

0207 

RDTX  -  E< 1 )**EX 

ISN 

0208 

THT!!  -  EX*THT! 1 > 

ISN 

0209 

US  «  RDTX*DCOS( THTX  ) 

ISN 

0210 

VS  -  RDTX*DSIN!THTX) 

ISN 

0211 

WRITE (6, 325  )  XS  ,  VS , US ,  VS , E ( 1 >,THT< 1  ), RDTX, THTX 

ISN 

0212 

232 

FORMAT! IS, 1P8E14.5  > 

ISN 

0213 

324 

FORMAT! '  LE  '.1P8E14.5) 

ISN 

0214 

325 

FORMAT!'  TE  '.1P8E14.5) 

ISN 

0215 

KOLM  ■  KJLE  -  1 

ISN 

0216 

DO  65  KJ  ■  2,  KOLM 

ISN 

0217 

K  -  KOLE  -  KJ 

ISN 

0218 

XP  -  RDP!K)*DCOS!THP!K)> 

ISN 

0219 

YP  -  RDPI K >*D3INITHP( K > > 

ISN 

0220 

UP  -  RDPX!K>*DCOS!THPX<K>  > 

ISN 

0221 

VP  -  RDPX!K)*DSIN!THPX!K) ) 

ISN 

0222 

WRITE!6,232)K0,XP.YP,UP,VP,RDP!K>,THP(K>,RDPX!K),THPX(K 

ISN 

0223 

65 

CONTINUE 

ISN 

0224 

XS  -  E<  KJLE  )*DCOS<  THT! KOLE ) ) 

ISN 

0225 

YS  -  E!KJLE  )*OSINITHT!KOLE ) > 

ISN 

0226 

RDL X  -  E<  KJLE  >**EX 

ISN 

0227 

THL X  -  EX*THT<  KOLE ) 

ISN 

0228 

US  -  R0LX*0C0S!THLX) 

ISN 

0229 

VS  -  RDLX*OSIN(THLX> 

ISN 

0230 

WRITE ! 6 , 324 )XS , YS , US , VS ,E( KOLE ) ,THT! KJLE  >,RDLX,THLX 

ISN 

0231 

KJMXM  -  KJMX  -  1 

ISN 

0232 

KJLP  -  KOLE  +  1 

ISN 

0233 

00  31  KJ  -  KJLP , KJMXM 

ISN 

0234 

K  •  KJ  -  KJLE 

ISN 

0235 

XS  -  RDS!K)*OCOS!THS!K>  > 

ISN 

0236 

YS  »  RDS!K)*DSIN!THS!K) ) 

ISN 

0237 

US  «  RDSX! K  )*DCOS<  THSX! K ) ) 

ISN 

0238 

VS  -  RDSX!K)*OSIN!THSX!K)> 

ISN 

0239 

WRITE! 6,232)  KJ , XS , YS , US , VS , RDS! K > ,THSI K ) , RDSX! K ) , THSX!  1 

ISN 

0240 

31 

CONTINUE 

ISN 

0241 

XP  -  E!KOMX)*DCOS!THT!KOMX>) 

ISN 

0242 

YP  -  E!KOMX)*DSIN! THT ! KJMX ) ) 

ISN 

0243 

RDTX  -  E! KJMX  )**EX 

ISN 

0244 

THTX  »  EX*THT( KJMX ) 

ISN 

0245 

UP  -  RDTX*DCOS! THTX ) 

ISN 

0246 

VP  -  ROTX*DSIN!THTX) 

ISN 

0247 

WRITE!6,325)XP,YP,UP,VP,E(K0MX), THT! KJMX), RDTX, THTX 

ISN 

0248 

WRITE! 6,242 >ZPLUS , ZMINUS 

ISN 

0249 

242 

FORMAT! //10X, 'POINTS  AT  INFINITY  ARE  LOCATED  IN  THE  CAP 

'  PLANE  AT: ' , 

//  15X , 'PLUS: ' , 1P2E15.4, ' 


MINUS: * ,2E15.4,//> 


DETERMINATION  OF  ZC  SUCH  AS  TO  MINIMIZE  THE  RATIO  RMAX/RMIN  IN  THE 
L.C.  OMEGA  PLANE 


<  ZMGA-ZMGB  >/( ZPLUS-ZMINUS ) 

<  ZMGA*ZMINUS-ZMGB*ZPLUS )/( ZPLUS-ZMINUS ) 

<  ZMGA*ZPLUS-ZMGB*ZMINUS )/( ZPLUS-ZMINUS  > 


ISN  0254 
ISN  0255 


ISN  0256 
ISN  0257 
ISN  0258 


ISN  0260 

ISN  0261 
ISN  0262 
ISN  0263 
ISN  0264 
ISN  0265 
ISN  0266 
ISN  0267 
ISN  0268 
ISM  0263 
ISN  0270 
ISN  0271 
ISN  0272 
ISN  0273 
ISN  0274 
ISM  0275 
ISN  0276 
ISN  0278 
ISN  0279 
ISN  0280 
ISN  0281 
ISN  0282 
ISN  0283 
ISN  0284 
ISN  0285 
ISN  0286 
ISN  0283 
ISN  0283 
ISN  0290 
ISN  0291 
ISN  0292 
ISN  0293 
ISN  0294 
ISN  0295 
ISN  0296 
ISN  0297 
ISN  0298 
ISN  0299 
ISN  0300 
ISN  0301 
ISN  0302 
ISM  0303 
ISN  0304 
ISN  0305 
ISN  0306 
ISN  0307 
ISN  0308 
ISN  0309 


WR ITE { 6 , SO  1  ) 

601  FORMAT <  ' 1 1  TER ' , 1 1 X . 1  ZD ‘ , 22X , 1  ZB * . 22X , 1 ZC 1 , 1 8X , 'ZOMSTR' ,  19X, ' ZNTRD 1 
1 / 1 3X , ' RMIN ' ,  8X , ' RMAX ' ,7X, 'RATIO'/'  ( ZA< KJ ) , KJ» 1 , KJMX > * > 

250  ITEP.-l 

RATIO-O.ODO 
I F  < IGOT.EQ. 1  )  GO  TO  60 
C 

C  FOR  A  FIRST  GUESS,  USE  ZO<  - 1 . 0  ,  + 1 . 0  ) 

C 

ZC«CCMPLX< -1 . 000,1.000) 

C 

60  ZB  -  <  ZMGA*ZPLUS-ZMGB*ZMINUS-ZC*( ZMGA-ZMGB  )  )/<ZPLUS-ZMINUS> 

ZD  -  <  ZMGB*ZPLUS-ZMGA*ZMINUS+( ZMGA-ZMGB  )/ZC)/(ZPLUS-ZMINUS) 

68  CONTINUE 

DO  75  K  -  1  ,KJS 
RS  -  RDSX( K  ) 

ARS  -  THSX<  K  > 

ZOMS(K)  -  DCMPLXf  RS*DCOS( ARS  >,RS*DSIN(ARS)> 

ZOMSUO  -  ( ZD-ZB*ZOMS (  K  >/ZC  )/ (  Z1 -ZOMS( K  )/ZC  > 

75  CONTINUE 

DO  85  K  -  l.KJP 
RP  -  RDPX(K) 

ARP  -  THPX(K> 

ZOMP(K)  =■  OCMPLXI  RP*DCOS<ARP  )  ,RP*DSIN(ARP  )  ) 

ZOMP <  K  >  -  ( ZD -ZB "ZOMP <  K  ) /  ZC  >/( Z1 -ZOMP  <  K  )  /  ZC  ) 

85  CONTINUE 

I F < ILE.EQ.l  )  GO  TO  II 
ZOMLE  -  ZB 
GO  TO  12 

11  RD  -  E(KOLE) 

TH  «  THT( KJLE  ) 

RS  -  RD"*EX 

TH  •  £ X*TH 

ZOMLE  -  DCMPLXI RS*OCOS( TH  )  ,RS*DSIN( TH  )  ) 

ZOMLE  -  (  ZO-ZB*ZOMLE /ZC  >/( Z1 -ZOMLE /ZC  ) 

12  I F < ITE .EQ.  1  )  GO  TO  32 
ZOMTE  -  ZD 

ZA( 1 )  -  ZOMTE 
GO  TO  33 

32  RD«E<  1  ) 

TH-THTI 1  ) 

RS-RD**EX 

TH-EX*TH 

ZOMTE-OCMPLX( RS*DCOS<  TH ) , RS*DS I N( TH  ) ) 

ZOMTE *<  ZD -ZB "ZOMTE /ZC )/ ( Z1 -ZOMTE /ZC  ) 

ZA( 1 ) -ZOMTE 

33  DO  76  KO  -  2.K0LM 
K  -  KJLE  -  KJ 

76  ZA<  KJ >  -  ZOMP(K) 

ZA(KJLE)  -  ZOMLE 
DO  77  K  •  l.KJS 
KJ  -  KJLE  ♦  K 

77  ZA(KJ)  -  ZOMS(K) 

ZA(KJMX)  -  ZA< 1  ) 

ZNTRD  -  DCMPLXI O.ODO.O.ODO) 

AREA  -  0.000 
RMIN-CDABSI ZA( 1 > ) 

RMAX-RMIN 


B-7 


L 


ISN 

0310 

ZMAX-ZA<  1  > 

ISN 

0311 

ZMIN-ZA!  1  ) 

ISN 

0312 

KOMXX  ■  1 

ISN 

0313 

KJMN-1 

ISN 

0314 

DO  78  KJ  •  2.KJMX 

ISN 

0315 

DAREA  -  DABS< OREAL (  ZA( KJ- 1 > >«DIMAG! ZA! KJ > >-DREAL( 
*  D I MAG !  ZA  <  K J - 1  )  )  ) / 2 . 0 D 0 

ISN 

0316 

ZBR  •  <ZA< KJ-1 )+ZA< KJ > >/3. ODO 

ISN 

0317 

ZNTRD  ■  ZNTRD  ♦  ZBR*OAREA 

ISN 

0318 

RABS-CDABS! ZA! KJ ) ) 

ISN 

0319 

IF <  RABS . GE . RMIN  )  GO  TO  79 

ISN 

0321 

RMI N-RABS 

ISN 

0322 

ZMIN-ZA(KJ  ) 

ISN 

0323 

KJMN-KJ 

ISN 

0324 

GO  TO  78 

ISN 

0325 

79 

IF(RABS.LE.RMAX)  GO  TO  78 

ISN 

0327 

RMAX-RABS 

ISN 

0323 

ZMAX»ZA<  KJ  > 

ISN 

0329 

KJMXX  •  KJ 

ISN 

0330 

73 

AREA  -  AREA  +  OAREA 

ISN 

0331 

RATIO-RMAX/RMIN 

ISN 

0332 

ZNTRD  «  ZNTRD/AREA 

ISN 

0333 

ZOMSTR-  ZC*(ZNTRD-ZD>/( ZNTRD -ZB  > 

ISN 

0334 

WR I TE  (  6 , 602 )  I TE R , ZD , ZB , ZC , ZOMSTR , ZNTRD , RMI N , RMAX 
1  (ZA(KJ).KJ-l.KJMX) 

ISN 

0335 

602 

FORMAT! /IS, IP 10E 12. 4/E  17.4, 2E 12.4/1 10E13.5) ) 

ISN 

0336 

IF( RATIO. LT.RTOL >  GO  TO  63 

ISN 

0338 

I F <  IGOT.EQ. 1 )  GO  TO  63 

ISN 

0340 

ITER  -  ITER  ♦  1 

ISN 

0341 

I F ( ITER. LE. 30)  GO  TO  62 

ISN 

0343 

WRITE ( 6 , 204 ) 

ISN 

0344 

204 

FORMAT! ///10X, ‘TOLERANCE  SPECIFIED  FOR  RMAX /RMIN  1 
*  1  30  ITERATIONS') 

ISN 

0345 

STOP 

ISN 

0346 

62 

CONTINUE 

ISN 

0347 

I F ! M . EQ . 2  >  GO  TO  66 

ISN 

0349 

ZDS- 1 . 1D0-ZMIN 

ISN 

0350 

KJ-KJMN 

ISN 

0351 

67 

RD-E! KJ  )**EX 

ISN 

0352 

TH-EX-THT! KJ ) 

ISN 

0353 

ZOM-DCMPLX! ROCOCOS !TH),RD-DS IN!  TH  > ) 

ISN 

0354 

ZC-! ZOM-! ZDS-ZG )+ZE  >/! ZDS+ZF -ZE *ZOM > 

ISN 

0355 

GO  TO  60 

ISN 

0356 

66 

M-l 

ISN 

0357 

ZDS-0 . 900*ZMAX 

ISN 

0358 

KJ-KJMXX 

ISN 

0359 

GO  TO  67 

ISN 

0360 

63 

IGOT  -  1 

ISN 

0361 

WR I TE ( 6 , 208  )  ZD , ZB , ZC , ZNTRD , ZOMSTR 

ISN 

0362 

208 

FORMAT! //10X,  ’CONSTANTS  FOR  MAPPING  FROM 

*  'Z  -  PLANE  TO  OMEGA  -  PLANE  ARE1. 

*  //20X , 1 A  -  1 , 1P2EZ0 . 5 , /20X , 1 B  -  1 , 1P2E20 . 5 , /20X , 1 C 

*  1P2E20.5.//20X, ‘ZNTRD  "  ’  ,  1P2E20 . 5 , /20X . 1 ZOMSTR  - 


1 .1P2E20.S) 


SET  UP  THE  ARRAYS  OF  THETA  AND  LN( R ) 


DO  41  KJ  -  l.KJMXM 

XI  KJ  )  -  DATAN2( D IMAG( ZA( KJ  )  )  , DREAL ! ZA( KJ  )  )  ) 


1SN  0363 
ISN  0364 


ISN  0365 
ISN  0366 
ISN  0367 


ISN  0363 
ISN  0369 
ISN  0370 
ISN  0371 
ISM  0373 
ISN  0374 
ISN  0375 

ISN  0376 
ISN  0377 
ISN  0378 
ISN  0379 
ISN  0381 
ISN  0383 
ISN  0384 
ISN  0386 
ISN  0387 
ISN  0386 

ISN  0389 
ISN  0390 
ISN  0391 
ISN  0392 
ISN  0393 


ISN  0394 

ISN  0395 
ISN  0396 
ISN  0397 
ISN  0398 
ISN  0399 
ISN  0400 
ISN  0401 
ISN  0402 
ISN  0403 
ISN  0404 
ISN  0405 

ISN  0406 


ISN  0407 
ISN  0408 


41  Y<KJ>  -  C0A8S < ZA( KJ ) ) 

X<  KJMX )  «  X< 1  )  +  TP  I 
Y <  KJMX  >  «  Y<  1  ) 

C 

C  NOW  ADJUST  THE  ARGUMENTS  OF  THE  THETA  ARRAY,  SO  AS  TO  BE  CONTINUOUS 
C  ACROSS  THE  BRANCH  CUT  < ALONG  THE  NEGATIVE  REAL  AXIS)  OF  THE 
C  0ATAN2  FUNCTION.  THIS  ADJUSTMENT  ASSUMES  THAT  THE  CONTOUR  IS 
C  TRAVERSED  IN  A  COUNTERCLOCKWISE  DIRECTION. 

C 

SR  •  O.ODO 
PO  -  X< 1 ) 

DO  410  KJ  -  2.KJMXM 

IF<DABS<X<KJ >-PO).GT. PI)  BR  -  I.ODO 

PO  «  X<  KJ  ) 

X<  KJ )  •  X<  KJ  )  ♦  BR*TP  I 
410  CONTINUE 
C 

RMIN  -  10.000 

RMAX  -  O.ODO 

DO  49  K  -  1 , KJMX 

I F <  Y<  K ) . LT . RMIN )  RMIN  »  Y( K  ) 

49  I F <  Y<  K  ) . GT . RMAX  )  RMAX  -  Y<  K ) 

WARSCH  -  DSQRT<  RMAX/RMIN )  -  I.ODO 
IF< WARSCH . LT . 0 . 3D0 )  OM  -  I.ODO 
WRITE! 6 , 202 ) 

WRITE <  6 , 243 ) 

243  FORMAT! 3X ,  'BLADE-SURFACE  IMAGE  IN  THE  OMEGA  PLANE:', 

*  /3X , ' KJ ' , 6X  ,  'REAL '  , IOX . ' I MAG ' , 12X , ' R ' , 9X , 'THETA' ,/ > 

DO  51  KJ  -  1 , KJMX 

WRITE! 6,270)  KJ , 2A< KJ  ) , Y! KJ  ) . X! KJ  ) 

51  CONTINUE 

00  43  KJ  •  1 , KJMX 
43  Y<  KJ  )  -  DLOG<  Y<  KJ  )  > 

C 

C  USE  FFT  TO  FIND  A<  N )  AND  B<  N  ) 

C 

C  FIRST  SET  UP  THE  CONSTANTS  FOR  THE  LN<  R ) ,  THETA  SPLINE  FIT 
C 

CALL  CISPLN<Y,X,E,F , KJMX ,1,128,1  ) 

C 

ZI  -  DCMPLX! 0.000,1. ODO  > 

N  -  64 

N2  -  128 

NP1  •  N  ♦  1 

NP2  *  N  +  2 

PBN  ■  PI/64. ODO 

NB2  »  N/2 

NB2P1  -  NB2  +  1 

ZNN  -  DCMPLX<DFLOAT<N), O.ODO) 

ZW2N  -  DCMPLX! DCOS! PBN  ),DSIN( PBN)  ) 

OMM  -  1.000  -  OM 
C 

CALL  THDGRK 
C 

C  FIND  ZETAA  AND  ZETAB 
C 

ZABST  -  ZMGA 
ZBBST  *  ZMGB 


B-9 


ISN 

0409 

RABST  «  1.000 

ISN 

0410 

RBBST  -  1 .ODO 

ISN 

0411 

DR  =  0.1D0 

ISN 

0412 

DTH  -  10.0D0*PI/180.0D0 

ISN 

0413 

THA  -  -6 . 0D0*DTH 

ISN 

0414 

DO  21  I  -  1,6 

ISN 

0415 

R  -  .4900  ♦  DR*DFLOAT< 1-1  ) 

ISN 

0416 

DO  22  0  -  1,13 

ISN 

0417 

TH  -  OTH*DFLOAT! J-l  )  +  THA 

ISN 

0418 

ZTAGS  ■  DCMPLX! R*DCOS!  TH>, R*DSIN(TH>) 

ISN 

0419 

CALL  OMETAIA.B.ZMG  , ZTAGS , ZTANSR , 65 , 1 . 0D-00 , l ) 

ISN 

0420 

RA  -  CDABS! ZMG-ZMGA  > 

ISN 

0421 

IF!RA.GT. RABST)  GO  TO  23 

ISN 

0423 

RABST  -  P.A 

ISN 

0424 

ZABST  -  ZTAGS 

ISN 

0425 

23 

ZTAGS  -  -ZTAGS 

ISN 

0426 

CALL  OMETA! A , B , ZMG  , ZTAGS , ZTANSR , 65 , 1 . 00-00 , 1  ) 

ISN 

0427 

RB  •  CDABS!  ZMG-ZMGB ) 

ISN 

0423 

IF( RB.GT. RBBST)  GO  TO  22 

ISN 

0430 

RBBST  -  RB 

ISN 

0431 

ZBBST  »  ZTAGS 

ISN 

0432 

22 

CONTINUE 

ISN 

0433 

21 

CONTINUE 

ISN 

0434 

ZTAGS  -  ZABST 

ISN 

0435 

M  -  0 

ISN 

0436 

CALL  OMETA< A, B.ZMGA, ZTAGS, ZTANSR, 65, 1 .00-05, M) 

ISN 

0437 

I F  <  M . EQ . 0  >  GO  TO  260 

ISN 

0439 

WRITE! 6,261 )  ZTAGS, RABST 

ISN 

0440 

261 

FORMAT! //5X, 'OMETA  FAILEO  TO  CONVERGE  FOR  ZETA 
»  /10X,‘ ZTAGS  -  *  ,  1P2E13.5. ’  RABST  »',E13.5> 

A:  '  , 

ISN 

0441 

STOP 

ISN 

0442 

260 

CONTINUE 

ISN 

0443 

ZETAA  -  ZTANSR 

ISN 

0444 

ZTAGS  »  ZBBST 

ISN 

0445 

M  -  0 

ISN 

0445 

CALL  OMETA!A, B.ZMGB, ZTAGS, ZTANSR, 65, 1 .00-05. M) 

ISN 

0447 

IF(M.EQ.O)  GO  TO  262 

ISN 

0449 

WRITE! 6,263)  ZTAGS. RBBST 

ISN 

0450 

263 

FORMAT! //SX, 'OMETA  FAILED  TO  CONVERGE  FOR  ZETA 

B:  '  , 

i 

•  /10X, 'ZTAGS  -  ' , 1 PZE l 3 . 5 , '  RBBST  «',E13.5> 

ISN 

0451 

STOP 

ISN 

0452 

262 

CONTINUE 

ISN 

0453 

r 

ZETAB  -  ZTANSR 

C  FIND  GAMMA,  ALPHA,  BETA,  AND  S  FOR  MAPPING  TO  ETA 

r 

-  PLANE 

ISN 

0454 

v> 

AP  »  CDABS! ZETAA  ♦  ZETAB) 

ISN 

0455 

AM  *  CDABS! ZETAA  -  ZETAB) 

ISN 

0456 

AB  -  CDABS! ZETAA*ZETAB> 

ISN 

0457 

CHY  -  ( 2 . ODO-AP*AP  *2 . ODO*AB*AB )/AM/AM 

ISN 

0458 

RT  •  DSQRT (CHY*CHY-l.ODO) 

ISN 

0459 

CA  «  OSQRT! OABS(CHY  +  RT)  ) 

ISN 

0460 

C8  -  OSQRT !0A8S!CHY-RT>) 

ISN 

0461 

SS  -  DMIN1 ! CA , CB  > 

ISN 

0462 

ZAL  -  (2.0D0*ZETAA*ZETAB*!SS*SS*( ZETAA-ZETAB  ) -ZETAA -ZETAB 

* 

’  /DCONJG! ZETAA ) )/< SS*SS«! ZETAA-ZETAB >*ZETAA+ZETAB-2. 000/ 

* 

'  DCONJG! ZETAA ) ) 

ISN 

0463 

ZBT  ■  ! 2 . ODO*ZETAA*ZETAB-ZAL*<  ZETAA+ZETAB  >  >/! ZETAA+ZETAB-; 

B-10 


I SN  0464 


*  *ZAL  ) 

ZGM  ■  SS*(  ZETAA-ZBT )/( ZETAA-ZAL  > 

C 

ISn  0465  415  WR ITE ! 6 , 245  )  ZABST 

ISN  0466  245  FORMAT! // 1  OX BEST  GUESS  FOR  ZETA  A  IS  ZABST  -  * .1P2E12.3) 

ISN  0467  WRI TE <  6 , 246  >  ZBBST 

ISN  0468  246  FORMAT! // 10X BE  ST  GUESS  FOR  ZETA  B  IS  ZBBST  -  '.1P2E12.3) 

ISN  0469  WR ITE ! 6 , 2 1 5  >  ZETAA.ZETAB 

ISN  0470  215  FORMAT!//  SX.'ZETAA  -  '.1P2E13.5,'  ZETAB  -  '.2E13.5) 

ISN  0471  WRITE ! 6 , 2 1 6 )  ZAL , ZBT , ZGM , SS 

ISN  0472  216  FORMAT!//  5X, 'ALPHA  -  MP2E11.3,'  BETA  -  \2E11.3, 

*  '  GAMMA  «  ',2E11.3,'  S  »  '  .2E11.3) 

C 

C  FINDING  THE  LOCATION  OF  BLADE-SURFACE  POINTS  IN  THE  ZETA  PLANE  ONLY 
C  INVOLVES  PHI! THETA ) ,  SINCE  R  -  1.  USE  SPLINE  INTERPOLATION 
C 

ISN  0473  PHI ( 129  )  -  PHI!  1  )  +  TPI 

ISN  0474  £! 129  >  •  E! 1  )  ♦  TPI 

ISN  0475  CALL  C I SPLN! PH  I , E , THT , F , 1 29 , 1 , 1 , 2  ) 

ISN  0476  CALL  C I SPLN! PH  I , E , X , F . 1 29 . 2 , KOMX , 2 > 

C 

C  WHEN  THE  THETA  VS.  PHI  CURVE  IS  VERY  STEEP,  IT  MAY  HAPPEN  THAT  THE 
C  SPLINE-FITTED  PHI  VS.  THETA  CURVE  IS  NOT  MONOTONIC:  CHECK  NOW  WHETHER 
C  THIS  HAS  HAPPENED,  AND  REPLACE  THE  SPLINE-FITTED  DATA  WITH  LINEAR 
C  INTERPOLATES  WHEREVER  IT  HAS. 

C 

ISN  0477  DO  810  K  ■  2.KJMX 

ISN  0478  IF!F!K).GE.F!K-1  )  )  GO  TO  810 

ISN  0480  DO  811  I  »  2,129 

ISN  0481  IF  <  E 1  I  ).GT.X!K> >  GO  TO  812 

ISN  0483  811  CONTINUE 

ISN  0484  I  -  129 

ISN  0485  812  IP  »  I 

ISN  0486  IM  -  1-1 

ISN  0487  CON  -< PHI! IP ) -PH I ! IM ) > / ! E ! I P >-E ! IM ) > 

ISN  0488  0  -  K 

ISN  0489  DO  815  JO  -  1,20 

ISN  0490  IF!X!J+00).GT.£!IP>)  GO  TO  816 

ISN  0492  815  CONTINUE 

ISN  0493  816  KP  -  0  ♦  JO  -  1 

ISN  0494  DO  817  OJ  *  1,20 

ISN  0495  IF!XI0-00>.LT.E!IM>)  GO  TO  818 

ISN  0497  817  CONTINUE 

ISN  0498  818  KM  ■  0  -  JJ+  1 

ISN  0499  DO  319  JJ  -  KM , KP 

ISN  0500  TMP  -  F ! 00  ) 

ISN  0501  F ! 0 J )  ■  PHI 1 IM )  +  CON*! X! 00 )-E I IM ) ) 

ISN  0502  WR ITE ! 6 , 8 1 3  )  OJ,F(JO>,TMP 

ISN  0503  813  FORMAT! 5X, 'NOTE:  SPLINE  FIT  REPLACED  BY  LINEAR  INTERPOLATION  IN', 

*  '  FINDING  PHI!', 13,')  -’,F10.5,'  OLD  PHI  «',F10.5> 

ISN  0504  819  CONTINUE 

ISN  0505  K  ■  KP 

ISN  0506  810  CONTINUE 

C 

ISN  0507  WRITE ( 6 , 2 1 7  > 

ISN  0508  217  FORMAT! //I OX, 'MAPPING  FROM  OMEGA  -  PLANE  TO  ZETA  -  PLANE:', 

*  //  4X, 'K' , 14X, 'OMEGA' ,26X, 'ZETA' ,//> 

ISN  0509  DO  54  K  -  l.KJMX 


B-ll 


ISN  0510 
ISN  0511 
ISN  0512 
ISN  0513 
ISN  0515 
ISN  0516 
ISN  0517 


ISN  0518 


ISN  0519 
ISN  0520 
ISN  0521 
ISN  0522 
ISN  0523 
ISN  0524 
ISN  0525 
ISN  0526 
ISN  0527 
ISN  0528 
ISN  0529 
ISN  0530 


ISN  0531 
ISN  0532 
ISN  0533 


ISN  0534 
ISN  0535 
ISN  0536 
ISN  0537 
ISN  0538 
ISN  0539 
ISN  0540 


ISN  0541 
ISN  0542 
ISN  0543 
ISN  0544 
ISN  0545 
ISN  0546 
ISN  0548 
ISN  0549 
ISN  0550 
ISN  0551 
ISN  0552 
ISN  0554 
ISN  0555 
ISN  0556 


R  -  DEXP ( Y! K  ) ) 

ZX  -  DCMPLX!R*DCOS(X(K)),R*DSIN!X!K))) 

ZYG  -  DCMPLX! DCOS<  F!K)),DSIN!F!K))) 

IF(K.EO.l.OR.K.EQ.KJMX)  ZYG  -  Z1 
ZCC!K)  -  ZYG 

54  WRI TE ( 6 • 2 1 8  >  K.ZX.ZYG 
218  FORMAT! 15, 1P4E15.5  > 

ZCC  NOW  CONTAINS  ZETA  ON  THE  BLADE  SURFACES 

WRITE! 6 , 202 ) 

NOW  DO  THE  MAPPING  FROM  THE  ETA  -  PLANE  TO  THE  KSI  TILDE  -  PLANE. 
WHERE  ETA/S  -  SNIKSI  TILDE) 

AK  •  SS*SS 
a kq  ■ 

AKP  -  DSQRT! 1 . ODO-AKQ ) 

AKM  -  AKP*AKP 

CALL  £LLPT!PI,AK,RL,1 ) 

TBK  -  2  .  ODO*RL 
CTR  -  -RL 
CAPK  -  RL 

CALL  ELLPT!PI,AKP,RL.l > 

CAPKPM  ■  RL 

WRITE! 6 ,222 )  AK , CAPK , AKP , CAPKPM 

222  FORMAT! //10X,  'COMPLETE  ELLIPTIC  INTEGRALS  OF  K  AND  K  PRIME', 

*  '  ARE  AS  FOLLOWS:', 

*  // 1  OX , ‘ K! ' , F 1 0 . 6  ,  ‘  >  -  '  ,  F10.6.5X, 'K! ' , F10.6, '  >  -  '.F10.6) 
IFIND-0 

DO  55  I  -  l.K-JMX 

ZA! I  )  -  ZGM*!ZCC! I >-ZAL )/!ZCC! I )-ZBT) 

ZA( I >  NOW  HOLDS  ETA 

ZTD  -  ZA!I)/SS 
TAU  -  DREAUZTD) 

DLT  -  01 MAG! ZTD ) 

TSQ  -  TAU«TAU 
DSQ  -  DLT*OLT 

ART  -  1.000  +  AKQ*! TSQ  +  DSQ) 

RT  -  DSQRT! ( 1 . ODO-AKQ*TSQ >*! 1 .  ODO-AKQ*TSQ >  ♦  AKQ*DSQ*< 2 . ODO* 

*  < 1 . ODO+AKQ*TSQ  ) 

*  +AKQ*DSQ )  ) 

BRQ  -  l.ODO  ♦  TSQ  ♦  DSQ 

BRT  -  DSQRT! I  1 .ODO-TSQ)*! 1 .ODO-TSQ)  ♦  DSQ*! DSQ+2 . ODO+2 . ODO*TSQ ) ) 
ALM  -  ! BRQ-BRT  >*<  ART-RT )/4 . ODO/AKQ/TSQ 
SGA  -  ! TSQ  +  DSQ  -ALM) 

SGA  -  SGA/!SGA+1 . ODO-ALM*AKQ*< TSQ+DSQ ) ) 

IF(TAU.EQ. 0.000)  GO  TO  403 
RTALM  -  DSQRT! ALM >*TAU/DABS< TAU) 

GO  TO  404 

403  RTALM  -  O.ODO 

404  SGN  -  l.ODO 

IF!DLT.LT. O.ODO)  SGN  -  -l.ODO 
RTSGA  -  SGN*DSQRT( SGA ) 

RTALM  ■  DARS IN! RTALM) 

RTSGA  -  OARSIN! RTSGA) 


B-12 


ISN  0557  CALL  ELLPT! RTALM , AK , RL , 0  ) 

ISN  0558  CALL  ELLPT< RTSGA , AKP ,AG , 0 > 

ISN  0559  IFIAG.GE.0.0D0)  GO  TO  57 

ISN  0551  AG  -  -AG 

ISN  0562  RL  — RL  -  TBK 

ISN  0563  57  ZAKI)  -  DCMPLX!  RL  , AG  ) 

ISN  0564  I F (  IFINO.EQ. 1  )  GO  TO  55 

ISN  0566  IFU.EQ.l)  GO  TO  55 

ISN  0568  IF!!RL-0REAL!ZA1< 1-1 >)).LT.O.ODO)  GO  TO  55 

ISN  0570  KEDGE  •  I 

ISN  0571  KGM  -  I  -  1 

ISN  0572  IFIND  •  1 

ISN  0573  55  CONTINUE 

C 

C  ZAKI)  NOW  HOLDS  KSI  HAT 
C 

ISN  0574  WRITE ( 6 , 2 1 9 ) 

ISN  0575  219  FORMAT! ////10X, ‘MAPPING  FROM  THE  ETA  -  PLANE  TO  THE  KSI  HAT  -  *, 

*•  PLANE’, 

*  ///  4X, ‘K‘ , 15X, 'ETA‘ ,25X, ‘KSI  HAT ',//> 

ISN  0576  DO  56  K  -  l.KJMX 

ISN  0577  WR ITE ( 6 , 2 1 8 )  K , ZA< K  ) , ZA1 < K > 

ISN  0578  56  CONTINUE 

C 

C  NOW  SET  UP  A  GRID  IN  THE  KSI-HAT  PLANE,  AND  MAP  IT  BACK 
C  TO  THE  Z  -  PLANE: 

C 

ISN  0579  WRITE! 6 , 202  ) 

ISN  0580  WRITE!  6 , 205  ) 

ISN  0581  205  FORMAT! /5X, ’MAPPING  OF  A  GRID  IN  THE  KSI-HAT  PLANE’, 

*  //3X , ’K  L ' ,  8X, ’KSI  HAT' , 

*  15X, ’ETA’ ,16X, ’ZETA’ ,16X, ’OMEGA’ , 

*  13X, ’Z  MAPPED’ ,14X, 'ZXY' ,//> 

ISN  0582  ZA2! I  )  -  ZTE 

ISN  0583  ZA2! KJMX  )  -  ZTE 

ISN  0584  ZA2! KJLE  )  -  ZLE 

ISN  0585  DO  91  KJ  -  2 ,  KJLM 

ISN  0586  K  -  KJLE  -  KJ 

ISN  0587  91  ZA2!KJ)  -  ZPIK) 

ISN  0588  DO  92  K  -  1 ,KJS 

ISN  0589  KJ  «  KJLE  ♦  K 

ISN  0590  92  ZA2! KJ  )  *  ZS!K) 

C 

C  THE  ZA2  ARRAY  NOW  HOLDS  THE  BLADE-SURFACE  COORDINATES,  IN  THE  ORDER 
C  KJ  -  l:  TE 

C  KJ  -  2, KJLM:  PRESSURE  SIDE,  FROM  TE  TO  LE 
C  KJ  -  KJLE:  LE 

C  KJ  -  KJLP , KJMXM:  SUCTION  SIDE,  FROM  LE  TO  TE 

C  KJ  -  KJMX:  TE  AGAIN 

C 

ISN  0591  KMXM1  -  KMX  -  1 

ISN  0592  LMXM1  -  LMX  -  1 

ISN  0593  HCPKPM  -  CAPKPM/2 , ODO 

ISN  0594  TWOCPK  -  2.0D0*CAPK 

ISN  0595  THCPK  -  3.0D0*CAPK 

ISN  0596  FCPK  -  4.000»CAPK 

ISN  0597  ZAG  •  ZAL *ZGM 

ISN  0598  EXINV  -  l.ODO/EX 
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ISN  0599 
ISM  0600 
ISN  0601 
ISN  0603 
ISN  0604 
ISN  0605 
ISN  0607 
ISN  0608 
ISN  0609 
ISN  0610 
ISN  0611 
ISN  0612 


ISN  0613 


ISN  0614 
ISN  0615 
ISN  0617 
ISN  0619 
ISN  0620 
ISN  0622 
ISN  0623 
ISN  0624 
ISN  0626 
ISN  0627 
ISN  0628 
ISN  0629 
ISN  0630 
ISN  0631 
ISN  0632 
ISN  0634 
ISN  0636 
ISN  0638 
ISN  0640 

ISN  0641 
ISN  0642 


ISN  0643 
ISN  0644 
ISN  0645 
ISN  0646 
ISN  0647 
ISN  0648 


ISN  0649 
ISN  0651 
ISN  0653 


DXIR  -  FCPK/DFL0AT1 KMXMI ) 

SHXTE  -  DREAL <ZA 11  1  >  > 

IF <  SHXTE . GT .CTR  >  SHXTE  -  SHXTE-FCPK 

SHK  «  -CAPKPM*CAPKPM/4.0D0/1THCPK  +SHXTE  ) 

SHKINV  -  1 . 0D0/SHK 

IF (  ISHEAR.EQ.O )  SHKINV  -  O.ODO 

DSHX  -  2.0D0/DFLOATIKMXM1 ) 

OSHY  »  1 .  ODO/DFLOAT ( LMXM1 ) 

ZEETE  -  ZA( 1 ) 

10  -  0 
K  -  1 
L  -  1 
C 

C  K  -  LOOP  STARTS  HERE 
C 

760  CONTINUE 
C 

C  USE  LINEAR  INTERPOLATION  TO  GIVE  A  FIRST  GUESS  AT  Z  ON  THE 
C  BLADE  SURFACE 
C 

XIR  -  SHXTE+DFL0AT1 K-l )*DXIR 
IF1XIR.LT. -THCPK  >  XIR  -  XIR  ♦  FCPK 
IF <  X IR . GT . CAPK  )  XIR  -  XIR  -  FCPK 
DO  371  I  •  KEDGE , KJMX 
IF(DREAL(ZA1( I  )  J.LT.XIR)  GO  TO  372 

371  CONTINUE 

DO  373  I  -  2.KGM 

IF<  DREAL (ZA1( I ) > . LT. XIR )  GO  TO  372 
373  CONTINUE 
I  -  KEDGE 

372  KA  -  I 

KB  »  I  -  l 

XIA  -  DREAL ( ZA1 { KA ) ) 

XIB  -  0REAL1ZA11KB)  ) 

IF<  KA . NE . KEDGE  >  GO  TO  380 
IF1XIA.LT. O.ODOIXIA  -  XIA  *  FCPK 

IF1XIB.LT. O.ODOIXIB  -  XIB  ♦  FCPK 

IF1XIR.LT. O.OOOIXIR  -  XIR  ♦  FCPK 

380  ZGSA  -  1 1XIR-XIA)*ZAZ1KB)*IXIB-XIR  )*ZA21 KA))/IXIB-XIA) 

C 

381  ZGS'ZGSA 

SHX  -  -1.000+DSHX*DFLOAT!K-1 > 

C 

CL-  LOOP  STARTS  HERE 
C 

770  CONTINUE 

SHY»DSHY*0FL0AT1  L-l  > 

XIM-HCPKPM*!  1 .ODO-SHY) 

XIR  -  1XIM-HCPKPM)**2 

XIR  -  SHXTE +  TWOCPK*! 1.0D0+SHX)+XIR*SHKINV 

ZX1  -  DCMPLXI X I R , X IM ) 

C 

C  BYPASS  IMAGE  CALCULATIONS  FOR  POINTS  THAT  FALL  ON  THE  IMAGES  OF 
C  PLUS  OR  MINUS  INFINITY  OR  THE  T.E. 

C 

1FI ISHEAR.EQ.O)  GO  TO  84 

1F1L .EQ. 1 . AND.K.EO. 1 )  GO  TO  73 

IF1L.EQ. t .ANO.K.EQ.KMX)  GO  TO  73 
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ISN 

0655 

84 

CONTINUE 

ISN 

0655 

IFiL.EQ.LMX.AND.K.EQ.l )  GO  TO  74 

ISN 

0658 

IF!L.EQ.LMX.AND.K.EO.KMX>  GO  TO  74 

ISN 

0660 

IF< IOE . EQ. 0  >  GO  TO  80 

ISN 

0662 

IF! L  .EQ. lMX .AND .K. EQ. KMXH )  GO  TO  81 

ISN 

0664 

C 

GO  TO  80 

« 

ISN 

0665 

C 

c 

73 

ZEETA  -  ZEETE 

ISN 

0666 

Z ETA  -  ( ZBT*ZE£TA-ZAG  )/<  ZEETA-7GM  ) 

ISN 

0667 

ZOMA  -  ZOMTE 

ISN 

0663 

ZFNL  -  DCMPIX(SHX.SHV) 

ISN 

0669 

ZXV  -  ZTE ‘ZGAMMA 

ISN 

0670 

GO  TO  710 

ISN 

0671 

74 

ZEETA  -  DCMPLXISS, 0.000) 

ISN 

0672 

ZETA  -  ZETAA 

ISN 

0673 

ZOMA  -  ZMGA 

ISN 

0674 

ZFNL  -  DCMPLX( SHX , SHV  ) 

ISN 

0675 

ZXV  -  2.0D0*ZA!LMX-1 >-ZA!LMX-2> 

ISN 

0676 

GO  TO  710 

ISN 

0677 

81 

ZEETA  -  DCMPLX< -SS.O.ODO) 

ISN 

0678 

ZETA  -  ZETAB 

ISN 

0679 

ZOMA  -  ZMGB 

ISN 

0680 

ZXV  -  2 . 0D0*ZA( LMX- 1 >-ZA(LMX-2) 

ISN 

0681 

ZFNL  -  DCMPLXCSHX.SHV) 

ISN 

0682 

c 

GO  TO  710 

ISN 

0683 

L 

c 

80 

CONTINUE 

ISN 

0684 

CALL  JCELFN!XIR,XIM,AKQ,AKM,RLS,AGS,l) 

ISN 

0685 

ZEETA  -  SS*DCMPLX(RLS,AGS) 

ISN 

0686 

ZETA  -  <  ZBT*ZEETA-ZAG  )/! ZEETA-ZGM ) 

ISN 

0687 

CALL  OMETA( A , B , ZOMA , ZETA , ZTANSR , 65,1. 0D-00 , 1 ) 

C 

c 

NOW  DO  THE  NEVTON-RAPHSON  ITERATION. TO  FIND  Z,  GIVEN  ZBOMK 

ISN 

0688 

c 

ZBOM-ZC*!  ZOMA-ZD )/<  ZOMA-ZB ) 

ISN 

0689 

RAD  -  CDABS(ZBOM) 

ISN 

0690 

ARG  -  DATAN2< D I MAGI ZBOM ) , DREAL( ZBOM ) ) 

ISN 

0691 

RADO  -  RAD**EXI NV 

ISN 

0692 

ARGO  -  EXINV*ARG 

ISN 

0693 

ZBOMK  -  DCMPLX! RADO*DCOS!ARCO), RADO*DSIN! ARGO ) ) 

ISN 

0694 

IT«  l 

ISN 

0695 

90 

Z  -  ZGS 

ISN 

0696 

ZZT  -  ( Z-ZT )/ZDN 

ISN 

0697 

ZT1  -  ZPI*ZZT 

ISN 

0698 

ZT2  -  CDSIN< ZT1  ) 

ISN 

0699 

ZZN  -  ( Z-ZN  )/ZDN 

ISN 

0700 

ZT3  ■  ZPI»ZZN 

ISN 

0701 

ZT4  -  CDS  INC  ZT3 ) 

ISN 

0702 

ZF-ZT2/ZT4-ZS0MK 

ISN 

0703 

87 

ZFPM  -  ZP I *CDS IN( ZTN ) /Z0N/ZT4/ZT4 

ISN 

0704 

ZNEW  -  Z  -  ZF/ZFPM 

ISN 

0705 

IF( KNR.EQ.O )  GO  TO  83 

ISN 

0707 

IF <  KOUNT . EQ . 0 )  WRITE!  6.211 > 

ISN 

0709 

211 

FORMAT! //5X, ‘DIAGNOSTIC  OUTPUT  OF  THE  NEWTON-RAPHSON  ITERATIONS 

f 


ISN  071  C 
ISN  0711 
ISN  0712 
ISN  0713 
ISN  0715 
1SU  0716 
ISN  0717 
ISN  0718 
ISN  0720 
ISN  0722 
ISN  072* 
ISN  0723 
ISN  0726 
ISN  0728 

ISN  0729 
ISN  0730 


ISN  0731 
ISN  0733 
ISN  0735 
ISN  0737 
ISN  0739 
ISN  0740 
ISN  07*2 
ISN  0743 
ISN  07*4 
ISN  0746 
ISN  0747 
ISN  0743 
ISN  0750 
ISN  0751 
ISN  0753 
ISN  0755 
ISN  0756 
ISN  0753 
ISN  0759 
ISN  0760 
ISN  0762 
ISN  0763 
ISN  0764 
ISN  0766 
ISN  0767 
ISN  0769 
ISN  0770 
ISN  0771 
ISN  0773 
ISN  0774 
ISN  0775 
ISN  0777 
ISN  0778 
ISN  078) 
ISN  0781 
ISN  0782 
ISN  0784 


*  /5X. 'VARIABLES  PRINTED  ARE  IT , ZBOMK , Z, ZN,EW, 2F , ZFPM ' , // ) 

WRITE ( 6,212  )  IT, ZBOMK. Z.ZNEV.ZF ,ZFPM 
212  FORMAT! 1 5 . 1 P 1 OE 12 . 3 > 

KOUNT  -  KOUNT  ♦  1 
IF(KOUNT.CT.INR)  STOP 
83  CONTINUE 
IT  -  IT  ♦  1 
ZGS  •  ZNEW 

I F ( CDABS ( ZNEW-Z  J.LT.1.0D-06)  GO  TO  72 
I F  <  CDABS  (  ZGS  > . GT  .  S I ZE  )  ZGS  -  ZGSA 
IF< IT.LE.50)  GO  TO  90 
82  ZGS  •  ZGSA 
ZNEW  -  ZERO 
I F  <  L . EQ . 1  )  ZNEW-ZGSA 
72  CONTINUE 
C 

ZFNL  -  DCMPLX! SHX , SHY  ) 

ZXV-ZGAMMA*ZNEW 

C 

C  CHECK  POINTS  ON  THE  PERIODIC  BOUNDARY,  NEAR  K-l,  K-KMX,  AND  K»KMX/2 
C 

I F ( ISHEAR.EQ.O  >  GO  TO  710 
IF(L.NE.LMX)  GO  TO  710 
IF(K.EQ.l.OR.K.EQ.KMX)  GO  TO  710 
IF( K . GT. 4  )  GO  TO  711 

714  DY  -  DIMAGIZXY-ZT) 

IF!DY.GE. 0.000)  GO  TO  717 
ZXY  -  ZXY  ♦  ZI *SG 

GO  TO  714 

717  IF(DY.LE.SG)  GO  TO  710 
ZXY  -  ZXY  -  ZI*SG 

GO  TO  714 

711  IF(K.GE.KAA)  GO  TO  712 
GO  TO  710 

712  IF <  K.GT. KMXH  )  GO  TO  713 
IF(K.EQ.KMXH.AND.IOE.EQ.l  )  GO  TO  710 

718  DY  ■  DIMAG(ZXY-ZN) 

IF ( DY .GE . 0 . ODO  )  GO  TO  719 
ZXY  »  ZXY  ♦  ZI  *SG 
GO  TO  718 

719  IF<  DY. LE . SG  >  GO  TO  710 
ZXY  -  ZXY  -  ZI *SG 

GO  TO  718 

713  IF(K.GT.KBB)  GO  TO  715 

720  DY  -  D IM4G! ZN-ZXY ) 

IF(DY.GE.O.ODO)  GO  TO  721 
ZXY  -  ZXY  -  ZI*SG 

GO  TO  720 

721  IF(DY.LE.SG)  GO  TO  710 
ZXY  -  ZXY  ♦  ZI *SG 

GO  TO  720 

715  IF<  K . LT . KMXM4  )  GO  TO  710 

722  DY  «  DIMAG! ZT-ZXY  ) 

IF! DY .GE . 0 . ODO  )  GO  TO  723 
ZXY  -  ZXY  -  ZI*SG 
GO  TO  722 

723  IF! DY .LE . SG  )  GO  TO  710 
ZXY  «  ZXY  -  ZI*SG 
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t 


I 


ISN 

0785 

C 

GO  TO  722 

ISN 

0786 

u 

c 

710 

WRITE <6, 224)  K.L.ZXI , ZEETA . ZETA , ZOMA , ZFNL , ZX V 

ISN 

0787 

ZAC l >  -  ZXY 

ISN 

0783 

XXSL.L)  -  DREAL(ZXY) 

ISN 

0789 

YYCK.L)  -  DIMAG(ZXY) 

ISN 

0790 

71 

CONTINUE 

ISN 

0791 

IFCIO.EQ.l)  GO  TO  790 

ISN 

0793 

L  •  L  +  1 

ISN 

0794 

IF(L.LE.LMX)  GO  TO  770 

ISN 

0795 

WRITEC  6 , 202 ) 

ISN 

0797 

I F <  PNCHZA  )  WR I TE ( 7 , 21 0 )  < ZAC L  ) , L-l , LMX ) 

ISN 

0799 

K  «  K  +  1 

ISN 

0800 

L  ■  1 

ISN 

0801 

IFCK.LE.KMX)  GO  TO  760 

ISN 

0803 

r* 

70 

CONTINUE 

V. 

c 

NOW  LOOK  FOR  CASES  WHERE  ZXY  -  ZERO,  AND  TRY  AGAIN,  USING 

c 

c 

INTERPOLATION  FROM  ALL  NEIGHBORING  POINTS 

ISN 

0804 

V# 

WRITE  <6,214) 

ISN 

0805 

214 

FORMAT! //5X, 'SECOND  ATTEMPT  TO  FIND  NON-CONVERGENT  CAS 

ISN 

0806 

WRITE <  6 , 205  > 

ISN 

0807 

DO  750  K  -  l.KMX 

ISN 

0808 

DO  755  L  -  2, LMX 

ISM 

0809 

IFCXXCK.D.NE.O.ODO)  GO  TO  755 

ISN 

0811 

IFC YYCK.L). NE. 0.000)  GO  TO  755 

ISN 

0813 

10  -  1 

ISM 

0814 

LI  *  L 

ISN 

0815 

K1  »  K 

ISN 

0816 

KM  •  K  -  1 

ISN 

0817 

KP  -  K  +  1 

ISN 

0818 

IFCKl.EQ.l )  KM  -  KMX  -  1 

ISN 

0820 

IF(Kl.EQ.KMX)  KP  -  2 

ISN 

0822 

LM  -  L  -  1 

ISN 

0823 

LP  -  L  ♦  l 

ISN 

0824 

Zll  ■  DCMPLXC  XXC  KM , LM ) , YYC  KM, LM ) ) 

ISN 

0825 

Z22  -  DCMPLXC XXCK1 ,LM) .YYCK1 ,LM) > 

ISN 

0826 

Z33  -  DCMPLXC  XXCKP.LM),YYCKP,LM)) 

ISN 

0827 

Z44  -  DCMPLXC  XXCKM.LI ),YYCKM,L1 ) > 

ISN 

0828 

Z55  -  DCMPLXC  XXC  KP , L  1 ),YYCKP,LI ) > 

ISN 

0829 

IFCU.EO.LMX)  GO  TO  756 

ISN 

0831 

Z66  -  DCMPLXC  XX(KM,LP),YY<KM,LP  >  > 

ISN 

0832 

Z77  -  DCMPLX<XX(KI,LP),YY(K1,LP)) 

ISN 

0833 

Z88  -  DCMPLXC  XX(KP,LP),YYCKP,LP)) 

ISN 

0834 

ZGSA  -  <Zll+Z22*Z33+Z44+Z55+Z66+Z77*Z88>/8.000 

ISN 

0833 

GO  TO  381 

ISN 

0836 

756 

ZGSA  •  (Z11*Z22+Z33+Z44+Z55)/5.0D0 

ISN 

0837 

GO  TO  381 

ISN 

0838 

790 

CONTINUE 

ISN 

0839 

755 

CONTINUE 

ISN 

0840 

750 

CONTINUE 

ISN 

0841 

STOP 

ISN 

0842 

100 

FORMAT! 8F10 . 4 ) 

ISN 

0843 

202 

FORMAT! Ill ) 

ISV 

0844 

210 

FORMAT! IP4E20. 13) 

ISN 

0845 

224 

FORMAT! 2I4,I2F10.5) 

ISN 

0846 

END 
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LEVEL  21.7  (  DEC  72  ) 


OS/360  FORTRAN  H 


COMPILER  OPTIONS  -  NAME-  MAI N , OPT-02 , L I NECNT-60 , S 1ZE -OOOOK , 

SOURCE, EBCDIC, NOLI  ST, NODECK, LOAD, MAP .NOEDIT. ID, XREF 
ISN  0002  SUBROUTINE  C I SPLN(  V  ,  X  ,  E  ,  F  ,  NP  ,  IRTN ,  NRTN  .  NPD  ) 

C 

C  THIS  SUBROUTINE  FITS  A  CUBIC  SPLINE  TO  A  FUNCTION  Y(X),  DEFINED  BY 
C  NP  PAIRS  OF  POINTS.  THE  SECOND  DERIVATIVE  OF  THE  FUNCTION 
C  IS  PERIODIC,  WITH  PERIOD  2-PI.  IF  NPD  -  I,  THE  FUNCTION 
C  ITSELF  IS  ALSO  PERIODIC,  WHILE  IF  NPD  -  2.  THE  FUNCTION  INCREASES 
C  BY  2 . -P I  EVERY  PERIOD:  Y(X+2.-PI>  -  Y  <  X  >  ♦  2. -PI.  SOLUTIONS  FOLLOW 
C  PAGES  9  -  15  OF 

C  THE  THEORY  OF  SPLINES  AND  THEIR  APPLICATIONS,  BY  0.  H.  AHLBERG , 

C  E.  N.  NILSON,  AND  0.  L.  WALSH,  ACADEMIC  PRESS,  1967 

C 

C 

C  NOTE  THAT  Y,X,£,AN0  F  ARE,  RESPECTIVELY,  ORDINATE,  ABSC I SSA , ABSC I SSA . 
C  AND  ORDINATE. 

C 

ISN  0003  IMPLICIT  REAL-81 A-H , O-Y > .COMPLE X- 1 6< Z > 

ISN  0004  DIMENSION  Y (  160), X(  160), E(  160), F( 1 60 ) , BDA(  160), EMI  160), H(  160) 

ISN  0005  DIMENSION  S< 1  SO ) , T< 1 60 ) , V< 1 60 ) , D< 1 60  > 

ISN  0006  DATA  TP  I /6 . 283 1 85307 1 79586/ 

C 

C  THIS  SECTION  (ENTERED  WHEN  IRTN  -  1)  USES  THE  NP  PAIRS  OF  INPUT 
C  COORDINATES  X  AND  Y  TO  FIND  THE  COEFFICIENTS  OF  THE  SPLINE  FIT. 

C 

C  THESE  COEFFICIENTS  -  HERE  CALLED  EM(KO)  -  ARE  THE  SECOND  DERIVATIVES 
C  OF  THE  FUNCTION. 

C 

C 


ISN 

0007 

I F < IRTN.EQ.2)  GO  TO  20 

ISN 

0003 

NPM  -  NP  -  1 

ISN 

0010 

N  -  NP  ♦  1 

ISN 

0011 

DO  1  KJ  -  2,NP 

ISN 

0012 

1  H(  K  1  )  -  X<  KO  )  -  X<  KO-1 > 

ISN 

0013 

H(N)  -  H ( 2  ) 

ISN 

0014 

DO  2  KJ  »  2 , NP 

ISN 

0015 

2  BDA(KJ>  »  H(K0  +  1  )/<H(KJ  )+H( KO+1 > ) 

ISN 

0016 

E<  1  )  -  O.ODO 

ISN 

0017 

F < 1  )  -  O.ODO 

ISN 

0015 

S(  1  )  -  1  .ODO 

ISN 

0019 

DO  3  KO  -  2, NPM 

ISN 

0020 

DN  -  2. ODO  +  (l.ODO  -  BDA< KO  )  >*E ( KO- 1 > 

ISN 

0021 

E ( KJ )  •  -BDA<  KJ  )/DN 

ISN 

0022 

0  (  K  J  >  - 

*  6.000*((Y<K0+1)-Y<K0>)/H(K0+1> 

*  -( Y<KO >-Y<  KJ-1 >)/H(KO))/(H(KO)+H(KO-l )> 

ISN 

0023 

S(  KJ  )  -  -(  1  .ODO-BOA(KJ)  )*S(KJ-1 1/DN 

ISN 

0024 

3  F ( KO  >  -  (D(KJ)-< 1.0D0-BDA(K0))*F(K0-l ) )/DN 

ISN 

0025 

Y(N)  -  Y<  2  1 

ISN 

0026 

I F ( NPD . EQ . 2  >  Y( N  )  -  YIN)  +  TPI 

ISN 

0028 

D(NP)  -  6.0D0-I (Y(N)-Y(NP)>/H(2> 

-  -<Y(NP)-Y(NPM))/H(NP>)/(H(NP)*H(2)> 

ISN 

0029 

NPMM  -  NP  -  2 

ISN 

0030 

T( NP  )  -  l.ODO 

ISN 

0031 

V( NP  )  ■  l.ODO 

ISN 

0032 

00  S  I  -  1 , NPMM 

ISN 

0033 

KO  -  NP  -  I 

ISN 

0034 

T(KJ)  -  E( KO )*T( KO  +  1  )  ♦  S(KO) 
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ISN  0035  V(  KJ  >  ■  E  (  KJ  )*V(  KJ  + 1  >  +  F(KJ) 

ISN  0035  6  CONTINUE 

ISN  0037  EM(NP)  »  (  D(  NP )-BOA<NP )*V<  2  )  -  ( 1 .ODO-BDA(NP  )  )*V( NPM)  >/ 

*  ( BDA( NP  )*T( 2  )  +  ( 1 .ODO-BDA(NP  >  )*T(NPM>  ♦  2.000) 

ISN  0038  DO  4  I  »  l.NPM 

ISN  0039  KJ  =  NP  -  I 

ISN  0040  4  EM( KJ  )  ■  E< KJ >*EM< KJ+1 >  +  FtKJ)  +  S(KO)*EM(NP) 

ISN  0041  SIZE  *  DA8S<  X<  NP  )  -  X< 1  )  > /  1 0  .  ODO 

ISN  0042  RETURN 

C 

C  THIS  SECTION  (ENTERED  WHEN  IRTN  -  2)  RETURNS  NRTN  INTERPOLATED 
C  VALUES  OF  THE  ORDINATE  F  AT  ASSIGNED  VALUES  OF  THE  ABSCISSA  E. 

C 

ISN  0043  20  KJ  =  2 

ISN  0044  DO  21  J  «  1 , NRTN 

ISN  0045  A  «  E ( J  ) 

ISM  0046  24  IF(A.LE.X(KJ  )  )  GO  TO  23 

ISN  0048  KJ  =  KJ  +  1 

ISN  0049  IF ( KJ . LE . NP  )  GO  TO  24 

ISN  0051  DF  =  A  -  X<  NP  > 

ISN  0052  IF  <  DF . GT . S IZE  )  WRITE(6,200  >  J ,  E ( J  > ,  NP  ,  X(  NP  > 

ISN  0054  200  FORMAT < / / 1  OX , ' WARN  I NG  -  ENTRY  IN  CISPLN  EXCEEDS  END  OF  BASE 

*  'ARRAY' ,/5X, *  E ( ' ,13, '  )  -  '.1PEI6.8,'  EXCEEDS  X(\I3,')  *  ', 

*  E16.8) 

ISN  0055  DF  »  A  -  X( 1 > 

ISN  0056  IF <  DF . LT . < -SIZE  )  )  WRITE ( 6 , 201  ) J , E ( J  ) , X( 1  ) 

ISN  0053  201  FORMAT(//10X, 'WARNING  -  ENTRY  IN  CISPLN  IS  LESS  THAN  THE  FIRST', 

*  '  BASE  POINT' , 

*  /5X, '£( 1  ,13, '  >  -  a,IPE16.8,'  IS  LESS  THAN  X( 1  )  »  '.E16.8) 

ISN  0053  KJ  =  NP 

ISN  0060  23  OXA  «  X(KJ)  -A 

I  ;,N  0061  DXB  *  A  -  X(  KJ-l  ) 

-!,f|  0062  CBA  -  DXA*DXA*0XA 

ISN  0063  CBB  -  DXB*DXB*DXB 

ISN  0064  F(J)  -  ( EM( KJ- 1 )*CBA+EM( KJ  >*CBB)/6.0D0/H(KJ  ) 

*  +  DXAMY(KJ-1  )-EM(KJ-l  )*H(KJ  >*H(K0)/6.0D0)/H(KJ  ) 

*  ♦  DXB*(Y<KJ )-EM(KJ )*H(KJ  >*H(KJ  >/6.0D0>/H(KJ  > 

ISN  0065  21  CONTINUE 

ISN  0065  RETURN 

ISN  0067  END 
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LEVEL  21.7  (  DEC  72  ) 


OS/360  FORTRAN  H 


COMPILER  OPTIONS  -  NAME-  MAIN .OPT-02 . L INECNT-60 , S IZE-0000K , 

SOURCE . EBCD IC , NOL I  ST . NODECK , LOAD , MAP . NOED IT , I D . XREF 
ISN  0002  SUBROUTINE  ELLPT( AR ,AK , ANS , KOMP  ) 

ISN  0003  IMPLICIT  REAL*8( A-H , 0-Y >. COMPLEX* 1 6< 2 ) 

ISN  0004  DIMENSION  SQI12) 

ISN  0005  DATA  K/l/ 

ISN  0006  DATA  P I /3 . 1 4 1 592653589793/ 

C 

C  THIS  SUBROUTINE  EVALUATES  THE  ELLIPTIC  INTEGRAL  OF  THE  FIRST  KIND, 
C  WITH  ARGUMENT  AR  (AN  ANGLE  IN  RADIANS),  AND  PARAMETER  AK  (A  REAL 
C  NUMBER). 

C  THIS  EVALUATION  USES  EQ.<U)  OF:  Y.  L.  LUKE,  ‘APPROXIMATIONS 
C  FOR  ELLIPTIC  INTEGRALS',  MATH.  COMP.,  VOL.  22  (JULY  1968),  PP  627- 
C  634,  WITH  N  ■  12. 

C  KOMP  -  0,1  FOR  THE  INCOMPLETE,  COMPLETE  INTEGRAL,  RESPECTIVELY. 

C 

ISN  0007  IF(K.GT.l)  GO  TO  11 

ISN  0009  K  -  2 

ISN  0010  TNP  ■  25.000 

ISN  0011  DO  10  M  -  1,12 

ISM  0012  THM  -  PI*FLOAT(M)/TNP 

ISN  0013  S  -  DSIN(THM) 

ISN  0014  10  SQ( M  >  -  S-S 

ISN  0013  11  AKK  -  AK*AK 

ISN  0016  SM  -  O.ODO 

ISN  0017  IF ( KOMP . EQ . 1  )  GO  TO  40 

ISN  0019  TN  -  DTAN(AR) 

ISN  0020  DO  20  M  -  1,12 

ISN  0021  SG  -  DSQRT( 1 . 0D0-AKK*SQ( M >  > 

ISN  0022  T  -  OATAN(SG-TN) 

ISN  0023  20  SM  *  SM  +  T/SG 

ISN  0024  ANS  -  (AR  +  2 . 0D0*SM )/TNP 

ISN  0025  RETURN 

ISN  0026  40  DO  41  M  -  1,12 

ISN  0027  SG  -  DSQRTI 1 . ODO-AKK*SO( M >  ) 

ISN  0028  41  SM  ■  SM  ♦  l.ODO/SG 

ISN  0029  ANS  -  PI*( 1 .000*2. ODO*SM >/2 . ODO/TNP 

ISN  0030  RETURN 

ISH  0031  END 
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LEVEL  21.7  <  DEC  72  ) 


OS/360  FORTRAN  H 


COMPILER  OPTIONS 


ISN  0002 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NAME-  MAIN.OPT-02.LINECNT-60.SIZE-0000K, 

SOURCE, EBCDIC, NOLI  ST, NODECK, LOAD, MAP, NOED IT, ID, XREF 
SUBROUTINE  JCELFNt RL , AG , AKQ, AKM , RLS , AGS , M ) 

WHEN  M.EQ. 1 , 

THIS  SUBROUTINE  RETURNS  THE  JACOBIAN  ELLIPTIC  SINE 
OF  A  COMPLEX  ARGUMENT: 

RLS  ♦  I *AGS  -  SNCRL  +  I*AG,K> 

USING  THE  ARITHMETIC  -  GEOMETRIC  MEAN  FORMULA  (SEE  P.571, 

OF  MATHEMATICAL  FUNCTIONS,  ED.  BY  M.  ABRAMOVITZ  AND  I.  A. 

U.S.  NATIONAL  BUREAU  OF  STANDARDS,  APPLIED  MATHEMATICS  SERIES,  55, 
JUNE  1964)  AND  THE  ADDITION  FORMULA  FOR  THE  SN  (SEE  EQUATION  125.01 
P.  24,  OF  HANDBOOK  OF  ELLIPTIC  INTEGRALS  FOR  ENGINEERS  AND 
PHYSICISTS,  BY  P.  F.  BYRD  AND  M.  D.  FRIEDMAN,  SPRINGER  VERLAG,  1954). 
AKQ  -  K**2,  AKM  -  1.  -  K**2 

WHEN  M.EQ. 2,  THE  QUANTITIES  RETURNED  ARE  THE  REAL  AND  IMAGINARY  PARTS 
OF  THE  PRODUCT  CN(  .  .  )*0N(  .  .  )  ,  WHICH  IS  THE  DERIVATIVE  OF  THE  SN(  . .  > 


HANDBOOK 

STEGUN, 


ISN 

0003 

IMPLICIT  REAL *8( A-H, 0-Y  >, COMPLEX* 16(Z> 

ISN 

0004 

DIMENSION  A(20),B(20),C(20),PH(20) 

ISN 

0005 

K  -  1 

ISN 

0006 

A(  1  )  -  1.000 

ISN 

0007 

B(  1  )  -  DSQRT(AKM) 

ISN 

0008 

C( 1 )  »  DSQRT ( AKQ ) 

ISN 

0009 

5 

DO  6  I  -  2,20 

isn 

0010 

All)  -  ( A ( 1-1 )+B( 1-1 >  >/2.0D0 

ISN 

0011 

B(  I  >  -DSQRT(A( 1-1 ) *B( 1-1 ) ) 

ISN 

0012 

C(  I  )  -  ( A < 1-1 )-B< 1-1 > J/2.0D0 

ISN 

0013 

IF(DABS(C< I ) ).LT. 1 .00-09)  GO  TO 

7 

ISN 

0015 

6 

CONTINUE 

ISM 

0016 

WRITE( 6 , 200 )  RL.AG 

ISN 

0017 

200 

FORMAT! ///10X, 'JCELFN  FAILED  TO 

CONVERGE 

ISN 

0018 

STOP 

ISN 

0019 

7 

NM  -  1-1 

ISN 

0020 

N  -  I 

ISN 

0021 

I F ( K . EQ . 2 )  GO  TO  20 

ISN 

0023 

PH(N)  -  A( N  >-RL*2**NM 

ISN 

0024 

IS 

DO  11  L  -  1 , NM 

ISN 

0025 

J  -  N  -  L 

ISN 

0026 

11 

PH(J)  -  ( PH( J+l  )+DARSIN( C( 0  +  1 )*1 

DSIN(PH( J 

ISN 

0027 

IF(K.EQ.Z)  GO  TO  40 

ISN 

0029 

SNK  -  DSIN( PH< 1 ) ) 

ISN 

0030 

CNK  -  DCOS(PH( 1 > ) 

ISN 

0031 

DNK  -  CNK/DCOS(PH( 2)-PH( 1 ) > 

ISN 

0032 

K  -  2 

ISN 

0033 

TMP  -8(1) 

ISN 

0034 

B(  1  )  -  C (  1  > 

ISN 

0035 

C(  1  )  -  TMP 

ISN 

0036 

GO  TO  5 

ISN 

0037 

20 

PH( N  )  -  A(N)*AG*2**NM 

ISN 

0033 

GO  TO  15 

ISN 

0039 

40 

SNP  -  OS I N( PH( 1  )  ) 

ISN 

0040 

CNP  -  OCOS(PH( 1 ) ) 

ISN 

0041 

DNP  -  CNP/DCOS( PH(2)-PH(1)) 

ISN 

0042 

DNM  -  1.000  -  SNP*SNP*DNK*DNK 

ISN 

0043 

IF( M. EQ. 2 )  GO  TO  50 

ISN 

0045 

RLS  ■  SNK*ONP/DNM 

FOR  Z  -  '  .1P2E15.4) 
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ISN 

0046 

AGS  -  CNK*DNK*SNP*CNP/DNM 

ISN 

0047 

RETURN 

ISN 

0043 

50  RLA  -  CNK*CNP 

ISN 

0049 

AGA  -  SNK*DNK*SNP*DNP 

ISN 

0050 

RLB  »  DNK*CNP *0NP 

isn 

0051 

AGB  -  AKG*SNK*CNK*SNP 

ISN 

0052 

RLS  -  <  RLA*RLB-AGA*AGB  l/DNM/DNM 

ISN 

0053 

AGS  -  < -RLA*AGB-RLB*AGA l/DNM/DNM 

ISN 

0054 

RETURN 

ISN 

0055 

END 
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T 


LEVEL  21.7  (  DEC  72  ) 


OS/360  FORTRAN  H 


ISN  0002 
ISN  0003 
ISN  0004 


COMPILER  OPTIONS  -  NAME-  MAIN , 0PT-02 , L I NECNT-60 , S IZE -OOOOK , 

SOURCE, EBCDIC, NOLI  ST, NODECK, LOAD. MAP. NOEDIT, ID, XREF 
SUBROUTINE  OMETA( A , B , ZMGA , ZTAGS , ZTANSR , N , E P S , M  ) 

IMPLICIT  REAL*8(A-H,0-V) , COMPLEX*! 6< Z ) 

DIMENSION  A(65),B(65),ZC(65) 


IF  M.EQ.O, 

THIS  SUBROUTINE  USES  NEWTON  -  RAPHSON  TO  FIND  ZETA( OMEGA  > :  ZMGA  IS  A 
KNOWN  VALUE  OF  OMEGA,  ZTAGS  IS  THE  INITIAL  GUESS  AT  ZETA,  ZTANSR  IS 
THE  SOLUTION.  AND  EPS  IS  THE  TOLERANCE  ON  THE  ANSWER. 


IF 

FOR 


M.EQ.l.,  THIS  SUBROUTINE  RETURNS  THE  VALUE  OF  OMEGA  (IN  ZMGA) 
A  GIVEN  VALUE  OF  ZETA  (IN  ZTAGS). 


IF  M.EQ.2,  THE  QUANTITY  RETURNED  (IN  ZTANSR)  IS  D  OMEGA/D  ZETA,  FOR 
GIVEN  VALUES  OF  OMEGA  (IN  ZMGA)  AND  ZETA  (IN  ZTAGS) 


ISN 

0005 

Z1  -  DCMPLXf 1.000,0 . ODO ) 

ISN 

0006 

ZETA  -  ZTAGS 

ISN 

0007 

NM  -  N  -  1 

ISN 

0003 

IT  -  1 

ISN 

0009 

5 

CONTINUE 

ISN 

0010 

ZSMA  -  DCMPLXf A( N  )  ,B( N  )  ) 

ISN 

0011 

IF(M.EQ. 1 >  GO  TO  22 

ISN 

0013 

ZSMB  -  ZSMA-DFLOAT(NM) 

ISN 

0014 

I F  <  M . EQ . 2  )  GO  TO  40 

ISN 

0016 

DO  10  J  -  1 , NM 

ISN 

0017 

ZSMA  -  ZETA-ZSMA  ♦  DCMPLXf A( N 

ISN 

0018 

ZSMB  -  ZETA-ZSMB  ♦  DCMPLXfAfN 

ISN 

0019 

10 

CONTINUE 

ISN 

0020 

ZEXP  -  CDEXP ( ZSMA  } 

ISN 

0021 

ZG  -  ZETA-ZEXP  -  ZMGA 

ISN 

0022 

ZDG  -  ZEXP *( Z1  +  ZSMB) 

ISN 

0023 

ZETOLD  -  ZETA 

ISN 

0024 

ZETA  -  ZETOLD  -  ZG/ZDG 

ISN 

0025 

IF  <  CDABSf ZETA -ZETOLD  > . LT. EPS  > 

ISN 

0027 

IT  -  IT  ♦  1 

ISN 

0028 

IF <  CDABSf  ZETA  > .GT . 1 . ODO  >  ZETA 

ISN 

0030 

I F ( IT.LE .50  )  GO  TO  5 

ISN 

0032 

M  -  5 

ISN 

0033 

RETURN 

ISM 

0034 

20 

ZTANSR  -  ZETA 

ISN 

0035 

RETURN 

ISN 

0036 

22 

DO  21  J  -  1 , NM 

ISN 

0037 

ZSMA  -  ZETA-ZSMA  +  DCMPLXf A(N 

ISN 

0038 

21 

CONTINUE 

ISM 

0039 

ZEXP  -  CDEXP ( ZSMA ) 

ISN 

0040 

ZMGA  -  ZETA-ZEXP 

ISN 

0041 

RETURN 

ISN 

0042 

40 

DO  41  0  -  1 , NM 

ISN 

0043 

41 

ZSMB  -  ZETA-ZSMB  +  DCMPLXf A(N 

ISN 

0044 

ZTANSR  -  ZMGA*( Z1 +ZSMB ) /ZTAGS 

ISN 

0045 

RETURN 

ISN 

0046 

END 

0.9D0*ZETA/CDABS(ZETA) 
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LEVEL 

21.7 

ISN 

1 

0002 

ISN 

0003 

ISN 

000  4 

ISN 

0005 

ISN 

0006 

ISN 

0007 

ISN 

0008 

ISN 

0009 

ISN 

0010 

ISN 

0011 

(  dec  72  )  OS/360  FORTRAN  H 

COMPILER  OPTIONS  -  NAME-  MAIN . OPT-02 . L INECNT-60 , S IZE-OOOOIC, 

SOURCE, EBCDIC, NOLI  ST, NODECK, LOAD, MAP, NOED IT. ID, XREF 
SUBROUTINE  SHAPE! C , H ,G , EX ,ZP , ZS , KOS , KOP  ) 

IMPLICIT  REAL  *8! A-H , 0- V ), COMPLEX- 1 6( Z ) 

DIMENSION  ZC( 80  )  ,ZP ( 80  ) 

DO  10  X  -  l.KOS 
10  READ! 5,100)  ZS!K> 

DO  20  K  -  1 , KJP 
20  READ! 5,100)  ZP!K) 

100  FORMAT! 8F 10.0) 

RETURN 

END 


'  LEVtC  2117't-DEC  ft  )  1  **  " *OS/380“_F OUTRAN  H 

COMPILER  OPTIONS  -  NAME-  MAI N , OPT-02 , L I NECNT-60 , S IZE-OOOOK , 

SOURCE, EBCDIC, NOLIST, NODECK, LOAD. MAP, NOEDIT. ID, XREF 
ISN  0002  SUBROUTINE  SHUFL! N.ZA.ZCC > 

C 

C  THIS  SUBROUTINE  TAKES  THE  N  COMPLEX  VALUES  IN  ARRAV  ZA,  WHICH  WERE 
C  COMPUTED  AND  STORED  IN  REVERSE  BINARY  ORDER  BY  FFT2  AND  "SHUFFLES* 
C  THEM  INTO  PROPER  ORDER  USING  ARRAY  ZCC  FOR  INTERMEDIATE  STORAGE. 

C  N  IS  ASSUMED  TO  HAVE  THE  FORM  2--M. 

C 


ISN 

0003 

IMPLICIT  REAL*8!A-H,0-Y) .COMPLEX* 16! Z ) 

ISN 

0004 

DIMENSION  ZA! 65 ) ,ZCC! 65 ) , IAL! 6 ) ,KR!  64 ) 

ISN 

0005 

DATA  KALL/O/ 

ISN 

0006 

DATA  IAL/6-0/ 

ISN 

0007 

IF  I KALL . EQ. 1 >  GO  TO  10 

IS'; 

0003 

KALL  -  1 

ISN 

•0010 

DO  341  JP  -  1 , N 

ISN 

0011 

0  -  JP  -  1 

ISN 

0012 

IAL! 6  )  -  0/32 

ISN 

0013 

0  -  0  -  32* IAL  <  6 ) 

ISN 

0014 

IAL! 5 )  -  0/16 

ISN 

0015 

0  -  J  -  1 6* IAL ( 5  ) 

ISN 

0016 

IAL ( 4 )  -  J/8 

ISN 

0017 

0  -  0  -  8*IAL! 4 ) 

ISN 

0018 

IAL ( 3  )  -  0/4 

ISN 

0019 

0  -  0  -  4*  IAL !  3  ) 

ISN 

0020 

IAL 12)  -  0/2 

ISN 

0021 

0  -  0  -  2*IAL! 2  > 

ISN 

0022 

IAL! 1 >  -  0 

ISN 

0023 

341 

1 

KR<  OP  )  - 

•  32* IAL! 1  )*16*IAL!2)  +  8*IAL(3)*4*IAL!4)*2*IAL!5)  +  IAL!6) 

ISN 

0024 

10 

DO  342  0  -  1 , N 

ISN 

0025 

342 

ZCC! 0  )  -  ZA! KR( 0  )♦  1  ) 

ISN 

0026 

DO  360  OP  -  1 , N 

ISN 

0027 

360 

ZA! OP )  -  ZCC! OP) 

ISN 

0028 

RETURN 

ISN 

0029 

END 

\ 
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LEVEL  21.7  (  DEC  72  > 


OS/360  FORTRAN  H 


ISN  0002 


ISN  0003 
ISN  0004 
ISN  0005 
ISN  0006 
ISN  0007 
ISN  0003 


ISN  0009 

ISN  0010 
ISN  0011 


ISN  0012 
ISN  0013 
ISN  0014 
ISN  0015 

ISN  0016 
ISN  0017 


ISN  0013 


COMPILER  OPTIONS  -  NAME*  MAI N .OPT-02 . L INECNT-60 , S IZE-OOOOK . 

SOURCE , EBCO IC . NOL I  ST, NODECK. LOAD , MAP , NOED IT, ID , XREF 
SUBROUTINE  THDGRK 


THIS  SUBROUTINE  MAPS  AN  OVAL  TO  A  UNIT  CIRCLE,  USING  A  VARIANT  OF 
THE  THECDORSEN-GARRICK  TRANSFORMATION  AND  FAST  FOURIER  TRANSFORM 
TECHNIQUES.  (SEE  REFERENCES  AT  BEGINNING  OF  MAIN  PROGRAM). 

IMPLICIT  R£AL*8!A-H,0-Y), C0MPLEX*1 6! Z ) 

COMMON/TGINTG/N , NP 1 , NP2 , N2 , NB2 , NB2P 1 , IP, IPMX, ITP , IWK, IMX.KJMX 
C0MM0N/TGCMPX/Z1 .ZW2N.ZI , ZNN , ZA , ZCC , ZA1 ,ZA2 
COMMON/TGDBLE/PBN , OM , OMM , ANGE RR , A , B , E , F , THT , PHI , X , Y 
DIMENSION  X( 160>,y( 160),E( 150),F( 150),THT( 160) 

DIMENSION  PHI <  1 30  )  ,A( 65  ) , B( 65 ) , ZCC( 65 ) , 

*  ZA( 165  )  ,ZA1( 165 > ,ZA2( 165 ) , 

*  IWK( 7 ) 

DIMENSION  ITP(IOO) 

301  DO  300  I  -  1.N2 

300  PHI!  I  >-PBN*DFLOAT< 1-1 ) 

THE  FOLLOWING  CALSPAN  LIBRARY  ROUTINE  PLACES  A  ZERO  IN  THE  LOCATIONS 
FROM  THE  FIRST  TO  THE  LAST  ARGUMENT  OF  THE  CALL.  THE  THIRD  ARGUMENT 
GIVES  THE  LENGTH  SPECIFICATION  OF  EACH  ENTRY. 


CALL  CLEAR(A< 1 >,A<65>,8> 

CALL  CLEAR( B( 1  ),B(65),B) 

CALL  CLEAR( E( 1 ),E( 150), 8) 

CALL  CLEAR! ZCC ( 1 ) , ZCC( 65 ) , 16  > 

WRITE! 6,401 > 

401  FORMAT! //I OX, ' PROGRESS  OF  PHI 


1  3X , 'IT' , 4X , ' OEMX ' , 4X , ' NO .  OF  THETA  REVERSALS') 
IT  -  1 

FIRST  GUESS  IS  THETA  -  THETA! TRAI L ING  EDGE)  -  PHI 


/  THETA  ITERATIONS  IS  AS  FOLLOWS:'// 


ISN  0019 
ISN  0020 
ISN  0021 
ISN  0022 


ISN  0023 
ISN  0024 
ISN  0025 
ISN  0026 
ISN  0027 
ISN  0028 
ISN  0029 
ISN  0030 
ISN  0031 
ISN  0032 
ISN  0033 
ISN  0034 
ISN  0035 


DO  315  K  -  1,N2 
315  E(  K  )  -  X!  1  )  ♦  PHUK) 

305  DEMX  -  0.000 
8! 1  )  •  0.000 
C 

C  USE  AJ  AND  BO  TO  GET  NEXT  APPROXIMATION  TO  THETA 
C 

ZCC(l)  -  DCMPLX! B( 1  ) , 0 . ODO  ) 

DO  302  OP  -  2.N 

302  ZCC(OP)  ■  DCMPLX ( B!OP)/2.ODO,-A!0P)/2.ODO) 
ZCC! NP 1  )  -  DCMPLX(0. 000,0. ODO) 

ZW  *  Z1/ZW2N 
DO  303  NP  -  1.NB2P1 

ZA 1 ( N P  >  *  ZCC! NP  )  +  DCONOG!ZCC(NP2-NP  ) ) 

ZW  »  ZW*ZW2N 

303  ZA2! NP  )  «  ZW* ! ZCC! NP >-DCONOG( ZCC! NP2-NP > )  ) 

DO  304  NP  -  1.NB2P1 

304  ZA! NP  )  ■  ZA1! NP  )  ♦  ZI*ZA2(NP) 

DO  309  NP  -  2.NB2 

309  ZA! NP2-NP  >  - 
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ISN  0036 
ISN  0037 
ISN  0038 
ISN  0030 
ISN  0040 
ISN  0041 
ISN  0042 
ISN  0043 
ISN  0044 
ISN  0046 
ISN  0047 
ISN  0048 
ISN  0049 
ISN  0050 
ISN  0051 
ISN  0053 
ISN  0054 
ISN  0055 
ISN  0057 
ISN  0058 
ISN  0059 

ISN  0060 
ISN  0061 
ISN  0062 
ISN  0063 
ISN  0064 


ISN  0066 


*  DCONJG! ZA1 ( NP  )  )  +  ZI*DCONJG< ZA2( NP  >  > 

CALL  FFT2< ZA , 6 , IWK ) 

CALL  SHUFL ( N , ZA , ZCC  ) 

B< 1  )  •  X<1)  -  PHI  (  1  )  -  OREAL! ZA(  I  )  ) 

DO  306  K  -  1,64 
TMP  -  E<  2*K-1 > 

E!2*K-1)  -  OR£AL(ZA(lO)  ♦  PHI!2«K-1)  ♦  B!  1  ) 

THT!2*K-1)  -  E!2*K-1> 

OEM  -  DABS!TMP-E!2*K-1>) 

IF ( DEM.GT . DEMX )  DEMX  -  DEM 
E!2*K-1)  -  0M*E!2*K-1)  ♦  OMM*TMP 
TMP  -  E!2*K) 

£{2"(C)  -  DIMAGI  ZA<  K  >  )  ♦  PHI  <  2»K  >  ♦  BID 

THTI 2*K )  -  E I 2*K ) 

DEM  -  DABS! TMP-E!2*K>) 

I F < DEM.GT. DEMX )  DEMX  -  DEM 
E! 2*K  )  -  0M*E!2*K)  ♦  OMM*TMP 
306  CONTINUE 

IF! IT.NE.ITPI 1P>)  GO  TO  840 
IP  *  IP  ♦  1 
WRI TE<  6 , 341  )  IT 
841  FORMAT! /3X. 

*  ‘THETA  BEFORE  AND  AFTER  RELAXATION  AT  IT  -  ’,14) 

WRITE! 6, 832)1 THTI I  ), 1-1 , 128) 

WRITE!  6, 832  HE!  I  ),  1-1,128) 

832  FORMAT! 1P10E13.5) 

840  CONTINUE 

IF! IT.GE.ITP! IPMX))  STOP 
C 

C  NOW  USE  THESE  THETAS  TO  GET  THE  NEXT  APPROXIMATION  TO  LN  R1K) 

C 

CALL  CISPLN!V,X,E,F, KJMX , 2 ,128,1) 

C 

C  NOW  FIND  THE  AJ  AND  00  COEFFICIENTS  CORRESPONDING  TO  THE  LN  RIK)  DATA 
C 


ISN 

0067 

DO  310  JP  -  1 , N 

ISN 

0063 

310 

ZA!  JP  )  -  DCMPLX!  F  ! 2*JP- 1  ) , F! 2*JP  )  ) 

ISN 

0069 

DO  307  NP  -  1 , N 

ISN 

0070 

307 

ZA! NP  )  -  DCONJG! ZA! NP  ) ) 

ISN 

0071 

CALL  FFT2!ZA,6, IWK) 

ISN 

0072 

DO  308  OP  ■  1 , N 

ISN 

0073 

308 

ZA! JP  )  -  DCONJG! ZA! JP  > )/ZNN 

ISN 

0074 

CALL  SHUFL ! N , ZA , ZCC ) 

ISN 

0075 

ZA! 65  )  -  ZA! 1  ) 

ISN 

0076 

DO  311  NP  >  l.NBZPl 

ISN 

0077 

ZA1 ! NP  )  -  ! DCONJG! ZA! NP2-NP  ) )  ♦  ZA!NP>)/2 

ISN 

0078 

311 

ZA2! NP  )  ■  ZI*IDC0NJC«ZA!NP2-NP>)  -  ZAINP) 

ISN 

0 079 

ZW  -  ZW2N 

ISN 

0080 

DO  312  NP  -  1.NB2P1 

ISN 

0081 

ZW  -  ZW/ZW2N 

ISN 

0082 

312 

ZCC! NP  )  -  !ZA1!NP 1  +  ZA2INP )*ZW)/2.0D0 

ISN 

0083 

ZW  *  Zl/ZW 

ISN 

0084 

DO  313  I  -  2.NB2 

ISN 

0085 

ZW  -  ZW-ZW2N 

ISN 

0086 

NP  *  NB2  ♦  I 

ISN 

0087 

ZBB  -  !  ZA1! NP2-NP )+ZA2! NP2-NP  >*ZW)/2 . ODO 

ISN 

0088 

313 

ZCC! NP  )  -  DCONJG! ZBB) 

ISN 

0089 

A!  1  )  ■  OREAL <  ZCC! 1 > ) 

ISN 

0090 

ZCC!NP1  )-0.5D0»DCONJG!ZAl( 1 >-ZA2< 1 >  > 

ISN 

0091 

A(  65  >  -  DREAL ( ZCC ( 65 ) ) 

ISN 

0092 

00  314  NP  »  2,N 

ISN 

0093 

A(NP)  -  2.000*DREAL!ZCC(NP ) ) 

ISN 

0094 

314 

B<  NP  >  -  2 . ODO*OIMAG ( ZCC ( NP ) > 

c 

c 

CHECK  FOR  CONVERGENCE 

ISN 

0095 

C 

I F < IT.EQ. 1 )  DEMX  -  1  .ODO 

ISN 

0097 

NRV  «  0 

ISN 

0098 

00  309  LL  •  2,128 

ISN 

0099 

809 

IF(E<LL  ).LE.E!LL-1  )>  NRV  »  NRV  *  1 

ISN 

0101 

WRITE! 6 , 400 )  IT, DEMX, NRV 

ISN 

0102 

400 

FORMAT < I5.1PE12.3,  I8> 

ISN 

0103 

IF ( DEMX , LT .ANGERR  >  GO  TO  316 

ISN 

0105 

IT  «  IT  +  1 

ISN 

0106 

I F ( IT.LE.IMX)  GO  TO  305 

ISN 

0108 

WRITE ( 6 , 2 1 3  ) 

ISN 

0109 

213 

FORMAT! //5X, ‘ITERATIONS  FOR  A  AND  B  DID  NOT  CONVERGE 

ISN 

0110 

c 

STOP 

ISN 

0111 

v» 

316 

WRITE! 6 , 203 )  IT, OM. DEMX 

ISN 

0112 

203 

FORMAT! //10X, ‘A  AND  B  ITERATIONS  CONVERGED  AT  IT  -  ' 

1 

•  USING  OM  -  '  ,  F6 . 3  ,  '  MAXIMUM  ANGULAR  ERROR 

»  1  PE  1 2 . 3 , ‘  RADIANS') 

ISN 

0113 

WRITE! 6 , 220 ) 

ISN 

0114 

220 

FORMAT!//  4X , 1 K ' , 1 1 X , ‘THETA’ ,15X, ‘PHI‘ ,//) 

ISN 

0115 

DO  69  K  -  1,128 

ISN 

OHS 

69 

WRITE ! 6 , 22 1 )  K , E! K  ) . PHI ! K ) 

ISN 

0117 

r 

221 

FORMAT! 15 , 1P2E20 . 6 ) 

ISN 

0118 

V 

RETURN 

ISN 

0119 

END 
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Appendix  C 

DICTIONARY  OF  VARIABLES 


FORTRAN  SYMBOL 

ALGEBRAIC 

EQUIVALENT 

DEFINITION,  USE,  COWENTS 

A(J) ,  B(J) 

fls  •  Bi 

Eq.  2-20,  2-21 

AK 

i  -  S‘ 

Eq.  2-23 

AKP 

Eq.  2-30 

AKQ 

V 

- 

ALM 

X 

Eq.  2-31 

CAPK 

K  Ck) 

Eq.  2-34 

CAPKPM 

K  Y£'> 

Eq.  2-34 

DLT 

6 

Eq.  2-28 

DXIM,  DXIR 

Eq.  2-38,  4-4 

EX 

0  r 

See  Figure  2 

EXINV 

2  ~  n 

G,  H 

&,  H 

See  Figure  1 

HCPKPM 

'N 

V/ 

* 

H<\j 

- 

OM 

7 r 

Relaxation  factor  used  in 
4>  ,  e  mapping 

PBN 

N 

PHI 

4 

- 

SGA 

<r 

Eq.  2-31 

SS 

S 

Eqs.  2-23 

TAU 

r 

C-l 

Eq.  2-28 

FORTRAN  SYMBOL 

ALGEBRAIC 

EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

TPI 

2  re 

- 

XIR,  XIM 

Ra  ( £ )  ,  'cJ'vvi  ( £ ) 

- 

ZA(N) 

ff\  cn) 

Eq.  A-9;  used  elsewhere 
for  temporary  storage 

ZB 

k> 

Eq.  2-9 

ZC 

/C 

Eq.  2-9 

ZD 

ou 

Eq.  2-9 

ZI 

V=T 

- 

ZN,  ZT,  ZLE,  ZTE 

ZN,?T 

See  Figure  1 

ZAL,  ZBT,  ZGM 

*,  /S  .  * 

Eq.  2-22 

ZA1(N),  ZA2 (N) , 

A,  (n.)  ,  (Kt)  , 

See,  for  example,  Eq. 

ZCC  (N) 

A-9.  Also  used  for 

C(At) 

temporary  storage 

ZBOM 

-/I 

- 

ZBOMK 

-nK 

ZEETA 

- 

ZETA 

S 

- 

ZNTRD 

centroid  of  the  a) -plane 

Z^MSTR 

St  -  plane  image  of  ZNTRO 

ZOMS(K),  ZOMP(K) 

- 

ZPLUS,  ZMINUS 

_n,+  ,  J2' 

Eq.  1-6 

ZXI 

< 

- 

ZXY 

%  +  iff 

C-2 


\ 


% 


Appendix  D 

LISTING  OF  METRIC  GENERATOR  PROGRAM 


LEVEL  21.7  (  DEC  72  )  OS/360  FORTRAN  H 

COMPILER  OPTIONS  -  NAME*  MAI N , OPT-02 . L I NECNT-60 . SIZE-0000K , 

SOURCE, EBCDIC, NOLI  ST. NODECK, LOAD, MAP, NOEDIT, I D.XREF 
C 

C  PROGRAM  CIMTVL  -  A  PROGRAM  TO  CALCULATE  THE  METRICS  OF  A  COORDINATE 
C  TRANSFORMATION  FOR  TURBOMACHINERY  CASCADES.  FOR  DOCUMENTATION  SEE: 

C 

C  0.  P.  NENNI  AND  W.  J.  SAE ,  EXPERIENCE  WITH  THE  DEVELOPMENT  OF 

C  AN  EULER  CODE  FOR  ROTOR  ROWS ,  ASME  PAPER  83-GT-36,  MARCH  1983 

C  V.  J.  RAE,  A  COMPUTER  PROGRAM  FOR  THE  IVES  TRANSFORMATION  IN 

C  TURBOMACHINERY  CASCADES,  AFOSR  TR-81-0154,  ADA096416,  NOV  1980 

C  W.  0.  RAE,  MODIFICATIONS  OF  THE  IVES-L IUTERMOZA  CONFORMAL - 

C  MAPPING  PROCEDURE  FOR  TURBOMACHINERY  CASCADES,  ASME  PAPER 

C  83  -GT- 116,  MARCH  1983 

C  W.J.  RAE,  REVISED  COMPUTER  PROGRAM  FOR  EVALUATING  THE  IVES 

C  TRANSFORMATION  IN  TURBOMACHINERY  CASCADES,  CALSPAN  CORPORATION 

C  REPORT  7177-A-l,  JULY  1983 

C 

C  THE  PROGRAM  READS  A  CARD  DECK  CONTAINING  THE  COORDINATES  X(K,D  , 

C  Y<  K.L  )  {  K  -  l.KMX  j  L  -  l.LMX,  AND  FINDS  BY  FINITE  DIFFERENCES 
C  THE  METRICS  OF  A  TRANSFORMATION  TO  A  RECTANGLE! KSI , ETA ) .  IN  THIS 
C  RECTANGLE,  THE  IMAGE  OF  THE  BLADE  SURFACE  LIES  ALONG  ONE  SIDE,  AND 
C  THE  IMAGES  OF  THE  POINTS  AT  INFINITY  LIE  AT  CORNERS  OF  THE  RECTANGLE 
C  (SEE  THE  REFERENCES  ABOVE).  THE  METRICS  ARE  WRITTEN  ON  TAPE  I. 

C 


ISN 

0002 

IMPLICIT  REAL*8(A-H,0~Z) 

ISti 

0003 

DIMENSION  Q<  4 1 0 , 6  ) 

ISN 

000* 

READ!  5,100)  KMX.LMX 

ISN 

0005 

100 

FORMAT! 2014) 

ISN 

0006 

READ! 5,101  )  SG 

ISN 

0007 

r 

101 

FORMAT! 8F 10. 4) 

V 

c 

r 

SG 

IS  THE  SLANT  GAP  BETWEEN  BLADES 

ISN 

0008 

DO  10  K  -  l.KMX 

ISN 

0009 

KM1LMX  -  ! K- 1  )*LMX 

ISN 

0010 

LKA  -  KM1LMX  ♦  1 

ISN 

0011 

LKB  -  KMiLMX  +  LMX 

ISN 

0012 

10 

READ! 5,200)! «Q!LK,S  >,QllK,6> > , LK-LKA , LKB > 

ISN 

0013 

c 

200 

FORMAT! 1P4E20. 13) 

c 

r 

Q(LK,5)  CONTAINS  X(K,L),  Q<LK,6>  CONTAINS  Y(l 

ISN 

0014 

V 

KM*  .KMX- 1 

ISN 

0015 

LM-LMX-1 

c 

r 

— 

SET  X  AND  Y  AT  K-KMX  ANO  L-l 

ISN 

0016 

v> 

LK-KM*LMX» 1 

ISN 

0017 

Ql LK.5  )*Q( 1 ,5  ) 

ISN 

0018 

c 

Q(LK,6)*Q( 1,6) 

V 

c 

r 

— 

FIND  X  AND  Y  AT  K-l  AND  L*LMX 

ISN 

0019 

Q(  LMX  ,5  )*2D0*'3<  LM,5  )-Q(  LMX-2,5  ) 

ISN 

0020 

Q!LMX,6)"ZD0*Q(LM,6)-Cl(LMX-2,6> 

C 

C  —  FIND  X  AND  Y  AT  K-KMX  ANO  L-LMX 
C 
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1 SN  0021 
ISN  0022 
ISN  0023 
ISN  0024 
ISN  0025 

ISN  0026 
ISN  0027 

ISN  0023 
ISN  0029 
ISN  0030 
ISN  0031 
ISN  0032 
ISN  0033 

ISN  0034 
ISN  0036 
ISN  0033 
ISN  0040 


ISN  0042 
ISN  00 43 
ISN  0044 
ISN  0045 
ISN  0046 

ISN  0047 
ISN  0049 


ISN  0051 
ISN  0052 
ISN  0053 
ISN  0054 
ISN  0055 
ISN  0056 
ISN  0057 


ISN  0058 
ISN  0059 
ISN  0060 
ISN  0061 
ISN  0062 
ISN  0063 


ISN  0064 
ISN  0065 
ISN  0066 
ISN  0067 
ISN  0068 


LK-KMX-LMX 

Q<LK,5)-2D0*Q(LK-l,5>-Q(LK-2,5) 

Q<LK,6>-2D0*Q<LK-1 , 5 >-G< LK-2 , 6 > 
TDELXI-4DO/DFLOAT(KM> 

TDELZT-2DO/DFLOAT(LM> 

C 

DO  520  K-l.KMX 
KM1LMX-IK-1 >*LMX 
C 

00  510  L-l.LMX 
LK-KM1LMX+L 
LP 1  *LK  +  1 
LM1-IK-1 
KP 1 -LK  +  LMX 
KM1-LK-LMX 
C 

IF< K. EO. 1 )  GO  TO  501 
IF(X.EQ.KMX)  GO  TO  501 
IFIL.EO.l >  GO  TO  504 
IF(L.EQ.LMX)  GO  TO  505 
C 

C  —  K-2  TO  KM  AND  L-2  TO  LM  < CENTERED  DIFFERENCES) 
C 

DXDXI-<Q<KP1,5>-Q<KM1,5>)/TDELXI 
DVDXI-<Q(KP1,6)-Q<KM1,6>)/TDELXI 
OXDZT-C  Q( LP1,5)-Q(LM1,5))/TDELZT 
DYDZT-I QI LP 1 , 6 )-Q<  LM1 , 6 ) J/TDELZT 
GO  TO  509 
C 

501  IFIL.EO.l >  GO  TO  510 
IFIL.EO.LMX)  GO  TO  510 
C 

C  —  K-l  OR  KMX  AND  L-2  TO  LM 
C 

LK2-LMX+L 
LKKM-IKM-l )*LMXH 
DXDXI-(Q(LK2,5)-Q(LKKM,5  >  )/TDELXI 
DYDXI-(Q<LK2,6)-Q( LKKM, 6  )  J/TDELXI 
0XD2T-I 0(LP1,5)-Q(LM1 ,5) J/TDELZT 
DYDZT-I 0(LP 1 , S )-Q< LM1 , 6 ) >/TDELZT 
GO  TO  509 
C 

C  -  l-l  AND  K-2  TO  KM 

C 

504  LK1-KM1LMX+1 
LK2-LK1-1 
LK3-LK2+1 

DXDZT-< -3D0*Q<  LK1 , 5 )+400*Q( LK2,5 )-Q<  LK3 , 5 ) )/TDELZT 
OYDZT-I -3D0*Q( LK1 , 6 )+4D0*Q( LK2,6)-Q(LK3,6> l/TDELZT 
GO  TO  506 
C 

C  -  L-LMX  ANO  K-2  TO  KM 

C 

505  LKS-I  KMX-K  )*LMX*LM 
LKB-KM1LMX+LM 

DXDZT-I Q( LKS, 5 )-Q( LKB , 5 ) )/TDELZT 
DQ-SG 

IFIK.GT.KMXH)  DQ--SG 


t 
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ISN  0070 

ISN  0071 
ISN  0072 

ISN  0073 


ISN  0074 
ISN  0073 
ISN  0076 
ISN  0077 
ISN  0078 
ISN  0079 


ISN  0080 


ISN  0087 
ISN  0088 


0YDZT»<Oi LKS,6)+DQ-QILKB,6) >/TDELZT 
C 

506  DXDXI*( Q( KP 1 , 5  >-Q( KM1 .5  )  I/TDELXI 
DYDXI ■( Q(  KP  1 ,6)-Q(KMl,6) >/TDELXI 
C 

509  RDMI-<0X0XI*0Y0ZT-DXDZT*DYDXI > 

C 

C  STORE  DERIVATIVES  IN  THE  ORDER  OKS  I /DX , DKSI /DY , DETA/DX , DETA/DY 
C 

CKLK.l  I-DYDZT/RDMl 
Q<LK,2>— DXDZT/RDMl 
Q<LK.3>«-DYDXI/RDM1 
Q( LK  ,  4  >•  DXDXI/RDM1 

510  CONTINUE 
520  CONTINUE 

C 

C  STORE  DERIVATIVES  ON  TAPE  1 
C 

LKMX>LMX*KMX 

C  WRITE! 1)  KMX,uMX.LKMX,<<Q<I,J>,I»l,LKMX>.J-l,4) 


C 

STOP 

END 
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